To restate the problem from Fabio Marroni, there are six plays:
These plays must be ordered without violating the following constraints:
“teatro” must take place after “affitto”, because in “affitto” we explain some things that will be needed for “teatro”. They don’t need to be consecutive, though.
Then a list of consecutive combinations that are forbidden
“affitto” cannot be followed by “tassista” nor by “teatro”
“autobus” cannot be followed by “ristorante”
“eredita” cannot be followed by “tassista” nor by “autobus”
“ristorante” cannot be followed by “affitto”, “tassista”, “teatro” nor by “autobus”
“tassista” cannot be followed by “affitto”, “teatro” nor by “eredita”
“teatro” cannot be followed by “autobus”, “affitto”, nor by “ristorante”
This is similar to the optimisation problem of scheduling subtasks. Basically the entire night's performance is one product. The product has six inputs (the individual scenes). Some inputs are composed of others and thus must follow one another.
I like to set these problems up in two main steps:
Then I will use the Rglpk package to solve the system.
Let's start by using \( x_i \) to denote the position of the \( i^{th} \) play from the list. For example \( x_1 = 3 \), denotes Affito is the third play of the night. Obviously \( 1 \leq x_i \leq 6 \).
The additional constraints can be expressed as:
| Constraint | Math | |
|---|---|---|
| “teatro” must take place after “affitto” | \( x_6 > x_1 \) | |
| “affitto” cannot be followed by “tassista” nor by “teatro” | \( x_1 + 1 \neq x_5,\; x_1 +1 \neq x_6 \) | |
| “autobus” cannot be followed by “ristorante” | \( x_2 +1 \neq x_4 \) | |
| “eredita” cannot be followed by “tassista” nor by “autobus” | \( x_3 +1 \neq x_5,\; x_3 +1 \neq x_2 \) | |
| “ristorante” cannot be followed by “affitto”, “tassista”, “teatro” nor by “autobus” | \( x_4+1 \neq x_1,\; x_4+1 \neq x_5,\; x_4+1 \neq x_6,\; x_4+1 \neq x_2 \) | |
| “tassista” cannot be followed by “affitto”, “teatro” nor by “eredita” | \( x_5+1 \neq x_1,\; x_5+1 \neq x_6,\; x_5+1 \neq x_3 \) | |
| “teatro” cannot be followed by “autobus”, “affitto”, nor by “ristorante” | \( x_6 +1 \neq x_2,\; x_6 +1 \neq x_1,\; x_6+1 \neq x_4 \) |
We will use \( y_{i,j} \in \{0,1\} \) to denote that play \( i \) is in position \( j \). So \( y_{2,4} =1 \) denotes that Autobus is the fourth play of the night.
Our first two new constraints ensure that we don't have two plays at the same time:
\[ \begin{align}
\sum_{i=1}^{6}y_{i,j} &= 1, \text{ for } 1 \leq j \leq 6 \\
\sum_{j=1}^{6}y_{i,j} &= 1, \text{ for } 1 \leq i \leq 6
\end{align}
\]
We can actually translate between the binary and discrete representations by observing:
\[ \sum_{j=1}^{6}jy_{i,j} = x_i, \text{ for } 1 \leq i \leq 6
\]
Thus the constraint table translates as follows:
| Constraint | Math | |
|---|---|---|
| “teatro” must take place after “affitto” | \( \sum_{j=1}^{6}jy_{6,j} > \sum_{j=1}^{6}jy_{1,j} \) | |
| “affitto” cannot be followed by “tassista” nor by “teatro” | \( \sum_{j=1}^{6}jy_{1,j} + 1 \neq \sum_{j=1}^{6}jy_{5,j},\; \sum_{j=1}^{6}jy_{1,j} +1 \neq \sum_{j=1}^{6}jy_{6,j} \) | |
| “autobus” cannot be followed by “ristorante” | \( \sum_{j=1}^{6}jy_{2,j} +1 \neq \sum_{j=1}^{6}jy_{4,j} \) | |
| “eredita” cannot be followed by “tassista” nor by “autobus” | \( \sum_{j=1}^{6}jy_{3,j} +1 \neq \sum_{j=1}^{6}jy_{5,j},\; \sum_{j=1}^{6}jy_{3,j} +1 \neq \sum_{j=1}^{6}jy_{2,j} \) | |
| “ristorante” cannot be followed by “affitto”, “tassista”, “teatro” nor by “autobus” | \( \sum_{j=1}^{6}jy_{4,j}+1 \neq x_1,\; \sum_{j=1}^{6}jy_{4,j}+1 \neq \sum_{j=1}^{6}jy_{5,j},\; \sum_{j=1}^{6}jy_{4,j}+1 \neq \sum_{j=1}^{6}jy_{6,j},\; \sum_{j=1}^{6}jy_{4,j}+1 \neq x_2 \) | |
| “tassista” cannot be followed by “affitto”, “teatro” nor by “eredita” | \( \sum_{j=1}^{6}jy_{5,j}+1 \neq x_1,\; \sum_{j=1}^{6}jy_{5,j}+1 \neq \sum_{j=1}^{6}jy_{6,j},\; \sum_{j=1}^{6}jy_{5,j}+1 \neq x_3 \) | |
| “teatro” cannot be followed by “autobus”, “affitto”, nor by “ristorante” | \( \sum_{j=1}^{6}jy_{6,j} +1 \neq x_2,\; \sum_{j=1}^{6}jy_{6,j} +1 \neq x_1,\; \sum_{j=1}^{6}jy_{6,j}+1 \neq \sum_{j=1}^{6}jy_{4,j} \) |
RglpkNote that we have 36 decision variables at this point. Eventually though you will see we need an additional variables for a total of 51 variables Our constraint matrix therefore will have 51 columns. Note that columes 1 through 6 denote \( y_{1,j},y_{2,j},...,y_{6,j} \), generally \( y_{i,j} \) is in column \( \left( j-1 \right)\cdot 6 + i \).
The first example is the constraint that Affitto must be in one and only one position.
\[ 1 y_{1,1} + 1 y_{1,2}+1 y_{1,3}+1 y_{1,4}+1 y_{1,5}+1 y_{1,6}+0 y_{2,1}+0 y_{2,2}+...+0 y_{6,6} = 1 \\ \]
The first row of our constraint table must be rearranged so that all decision variables are on the left hand side of the equation.
\[ \begin{align} -\sum_{j=1}^{6}jy_{1,j} + \sum_{j=1}^{6}jy_{6,j} &> 0 \\ -1y_{1,1} + -1y_{1,2} + -1y_{1,3} + -1y_{1,4} + -1y_{1,5} + -1y_{1,6} + 0y_{2,1} + 0y_{2,2} + ... + 0y_{5,6} + 1y_{6,1} + 1y_{6,2} + 1y_{6,3} + 1y_{6,4} + 1y_{6,5} + 1y_{6,6} &> 0 \end{align} \]
Note that constraints of the type “\( \neq \)” are not allowed in Rglpk. We will have to turn each “\( \neq \)” into two constraints, one “\( < \)” and one “\( > \)”.
It's a bit cumbersome to put the full constraint matrix into \( \LaTeX \). I'll jump straight into constructing the matrix in R.
It will take 12 constraints to ensure that one and only one play is in each slot:
exclusionMatrix<-matrix(nrow=12,ncol=51,data=0)
#Perhaps someone could help me do this without a for loop?
for (iRow in c(1:6))
{
#Each play is in one spot
exclusionMatrix[iRow,(iRow-1)*6+c(1:6)]<-1
#Each spot has one play
exclusionMatrix[iRow+6,(c(1:6)-1)*6+iRow]<-1
}
#All of these must sum to 1
exclusionRhs<-rep(1,times=12)
#All of these are strictly equal
exclusionDir<-rep("==",times=12)
The first row of our constraint table will take 1 row in the matrix.
tableRow1Matrix<-matrix(nrow=1,ncol=51,data=0)
#Set play 1 (Affito) to be negative
tableRow1Matrix[1,c(1:6)]<- -1
#Set play 6 (Teatro) to be positive
tableRow1Matrix[1,(6-1)*6+c(1:6)]<- 1
#The right hand side is 0
tableRow1Rhs<-c(0)
#This constraint is greater than
tableRow1Dir<-c(">")
The rest of the table contains 15 “\( \neq \)” constraints. Therefore we need 30 more rows in our matrix. Additionally we will actually need a few more decision variables. Let's look at the first “\( \neq \)” constraint:
\[ \begin{align} \sum_{j=1}^{6}jy_{1,j} + 1 &\neq \sum_{j=1}^{6}jy_{5,j} \\ \sum_{j=1}^{6}jy_{1,j} + -\sum_{j=1}^{6}jy_{5,j} &\neq -1\\ 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} &\neq -1 \\ \text{Which can be divided into} \\ \end{align} \]
\[ \begin{align} 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} &> -1 \\ \text{OR} \\ 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} &< -1 \\ \end{align} \]
We will introduce another set of binary variables denoted \( z_k \). If the “\( > \)” is binding then \( z_1 = 1 \). If the “\( < \)” is binding then \( z_1 =0 \). So we can turn our “OR” condition into an “AND” condition:
\[ \begin{align} 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} +10\left(1-z_1 \right) &> -1 \\ 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} +10 + -10z_1 &> -1 \\ \text{SO} \\ 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} + -10z_1 &> -11 \\ \text{AND} \\ 1y_{1,1} + 2y_{1,2} + ... + 6y_{1,6} + 0y_{2,1} +0y_{2,2} + ... + 0y_{4,6} + -1y_{5,1} + -2y_{5,2} + ... + -6y_{5,6} + 0y_{6,1} + ... + 0y_{6,6} -10z_1&< -1 \\ \end{align} \]
The choice of 10 above was arbitrary. We simply needed a number of sufficient magnitude to make the constraints non-binding. Now we will need \( z_k, 1 \leq k \leq 15 \). We can code these constraints as follows,
# These plays may not come directly before...
precede<-c(1,1, # These are
2, # broken out into the same
3,3, # rows as the table.
4,4,4,4,
5,5,5,
6,6,6)
# These plays may not come directly after.
follow<-c(5,6,
4,
5,2,
1,5,6,2,
1,6,3,
2,1,4)
notEqualMatrix<-matrix(nrow=30,ncol=51,data=0)
notEqualRhs<-rep(c(-11,-1) ,times=15)
notEqualDir<-rep(c(">","<"),times=15)
for (iRow in c(1:length(follow)))
{
#The "less than" condition
notEqualMatrix[2*iRow-1,(precede[iRow]-1)*6+c(1:6)]<-c(1:6)
notEqualMatrix[2*iRow-1,(follow[iRow]-1)*6+c(1:6)]<- -c(1:6)
notEqualMatrix[2*iRow-1,36+iRow]<- -10
#The "greater than" condition
notEqualMatrix[2*iRow,(precede[iRow]-1)*6+c(1:6)]<-c(1:6)
notEqualMatrix[2*iRow,(follow[iRow]-1)*6+c(1:6)]<- -c(1:6)
notEqualMatrix[2*iRow,36+iRow]<- -10
}
Now comes the creative part. Since we are using a linear programming model we must have an objective function, i.e. some linear function of the decision variables that we are trying to maximize or minimize. There is no innate objective function in Marroni's problem. I think an interesting objective would be to see if we can get the plays to be as close as possible to alphabetical order.
Note that if the plays are in alphabetical order then
\[ \begin{align} x_i &= i, \text{ for } 1 \leq i \leq 6 \\ \sum_{j=1}^{6}jy_{i,j} &= i, \text{ for } 1 \leq i \leq 6 \\ \end{align} \]
So let's set our constraint function to some sort of error term. The coefficient for \( y_{i,j} \) will be \( \left(i-j\right)^2 \). Note that we don't care about the \( z_k \) terms so they will have zero as their coefficients.
obj<-rep(0,times=51)
# (The ith play ) - (The jth position)
obj[c(1:36)] <- (((c(1:36) %/% 6) +1 ) - (c(1:36) %% 6))^2
Now let's put together all our constraints and solve this thing!
library(Rglpk)
## Loading required package: slam
## Using the GLPK callable library version 4.47
fullMatrix<-rbind(exclusionMatrix,
tableRow1Matrix,
notEqualMatrix)
fullRhs<-c(exclusionRhs,
tableRow1Rhs,
notEqualRhs)
fullDir<-c(exclusionDir,
tableRow1Dir,
notEqualDir)
#All our variables are binary
types<-rep("B",times=51)
#We want to minimize our objective function
max<-FALSE
Alphabet<-Rglpk_solve_LP(obj=obj,mat=fullMatrix,dir=fullDir,rhs=fullRhs,types=types,max=max,verbose=FALSE)
For comparison let's also try to figure out the order that brings us furthest from alphabetical order by maximizing the error scores.
max=TRUE
antiAlphabet<-Rglpk_solve_LP(obj=obj,mat=fullMatrix,dir=fullDir,rhs=fullRhs,types=types,max=max,verbose=FALSE)
Finally let's look at the orders
#We only need the first 36 variables
Alphabet$solution[1:36]*rep(c(1:6),times=6)
## [1] 0 0 0 0 0 6 1 0 0 0 0 0 0 2 0 0 0 0 0 0 3 0 0 0 0 0 0 4 0 0 0 0 0 0 5
## [36] 0
antiAlphabet$solution[1:36]*rep(c(1:6),times=6)
## [1] 0 0 0 0 5 0 0 0 0 4 0 0 0 0 3 0 0 0 0 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [36] 6
You should see for the alphabetical order that play 1 (affitto) is in position 6. So the order from first to last for the alphabetical optimization is:
Notice that they all just shifted up one position. Now for the antialphabetical order:
Thanks! I hope you enjoyed the RPub! Check out my blog at http://www.quantifyingQuality.com