Banner

Wednesday, January 3, 2018

Kotlin and Linear Programming Part I - Binary Programming

kotlin_for_operational_planning_and_optimization

In November 2017, I spoke at KotlinConf proposing Kotlin as a data science platform. When I was planning the code examples for the presentation, there were many topics to choose from. A common example would be a “classic” machine learning solution processing images, text, or speech with a neural network. I probably could have done something with TensorFlow which would undoubtedly have been cool.

But as fascinating as that area of data science is, I strongly felt the need to do something different. I wanted to present a topic that is highly relevant and for unknown reasons gets shadowed. I wanted to present something practical, accessible, and excitingly useful all at the same time, which would be especially appropriate for a Kotlin event. That is why I chose linear/integer programming and presented a driver scheduling example using Kotlin and ojAlgo.

Linear programming (also called linear optimization) is an applied field of mathematics often used in operations research and planning. It attempts to find an optimal solution to a planning problem when a set of business constraints exist. Some areas where linear programming is critical include:

  • Transportation networks (airlines, trains, ground, etc)
  • Staff scheduling (nurses, doctors, workers, teachers)
  • Resource planning (classroom scheduling, manufacturing, farming, supply chain management)
  • Product mixing (financial portfolios, ingredient mixing, etc)
  • Game “AI” (sudoku, chess, tic-tac-toe, etc)

More succinctly, when you hear a business problem that includes objectives that say “minimize”, “maximize”, or simply find a solution that’s “feasible”, linear programming should come right out of your tool kit. As a matter of fact, a case can be made that constraint solving is a form of machine learning. In many situations, it may even be a better approach than neural networks or other forms of AI.

Linear programming is a quiet underdog and workhorse, used heavily in operation-driven industries with little fanfare. Maybe I have my biases since I work for an airline, but it surprises me how little attention it gets. Even for consumer mobile apps, there is a wide spectrum of apps that have never been explored. For instance, how many mobile apps are there that will automatically schedule staff shifts? There are certainly big, expensive enterprise applications out there, but I have not seen a $5 Android app that does this.

Linear vs Integer vs Mixed Programming

“Linear” means continuous, which more or less means a variable can be be solved as a decimal value. If I am scheduling a worker for a given day, the optimal shift start can solve to 12, 12.5 hours, and even 12.5000012 hours. The number is allowed to float freely anywhere on a continuous number line, rather than using fixed whole numbers. Pure linear programming may be desirable for “mixing problems”, such as food ingredients or financial portfolio optimization.

For scheduling problems however, pure linear programming is not enough. Shifts often are modeled as discrete “slots” which can be represented with whole numbers. For instance, I may have a 12-hour day and I need to assign workers for that day. I could arbitrarily decide a “slot” is 1 hour, which would result in 12 slots. I could make a slot 4 hours, which would result in 3 slots or “shifts”. My model then describes the problem as “Slot1, Slot2, Slot3”… and so on. This is known as integer, or “discrete” optimization.

If you do not want to optimize with partial units, integer programming is handy for that situation as well. For instance, I may not want to manufacture 100.5002 units of a product. It must be a whole number like 100 or 101. Integer programming is also often used for binary/boolean logic by limiting a variable to 1 or 0, which is what I will focus on heavily in this blog post.

You can mix linear and integer programming, and this is known as “mixed integer programming”. This basically means some variables being solved are allowed to be continuous, while others may be integers. The driver problem I presented at KotlinConf was a mixed problem, as the driver shift start/ends were continuous, but to prevent overlap I used binary variables that could be 1 or 0.

But again in this blog post, I am going to focus heavily on binary optimization. I am hoping to do a series of blog posts on linear programming and optimization, but I think this will make a good first topic.

Getting Set Up

To keep concepts focused and simple, I am going to present an abstract problem. However this can be applied to real-life examples especially for scheduling applications. I am going to use Kotlin and ojAlgo, so be sure to configure both in your Maven/Gradle dependencies.

Before we dive in, also add these helper functions for ojAlgo to act as an improvised DSL. We will also instantiate an ExpressionsBasedModel and our helper functions will automatically name variables and functions before adding them to it.

 
import org.ojalgo.optimisation.ExpressionsBasedModel
import org.ojalgo.optimisation.Variable
import java.util.concurrent.atomic.AtomicInteger
 
// declare ojAlgo Model
val model = ExpressionsBasedModel()
 
// custom DSL for  expression inputs, eliminate naming and adding
val funcId = AtomicInteger(0)
val variableId = AtomicInteger(0)
fun variable() = Variable(variableId.incrementAndGet().toString().let { "Variable$it" }).apply(model::addVariable)
fun addExpression() = funcId.incrementAndGet().let { "Func$it"}.let { model.addExpression(it) }

An Abstract Example

Say we have five letters: A, B, C, D, and E. Let’s represent this as a Kotlin enumerable.

enum class Letter {
    A,
    B,
    C,
    D,
    E;
}

Now let’s create five numbers: ONE, TWO, THREE, FOUR, and FIVE. Please do not confuse this with the Java Number type. We are defining our own. The value property will hold its int counterpart.

enum class Number(val value: Int)  {
    ONE(1),
    TWO(2),
    THREE(3),
    FOUR(4),
    FIVE(5);
 
    override fun toString() = value.toString()
}

This pattern we are about to learn is useful and can be used to assign university classes, employee schedule shifts, etc. For example, Letter could represent an employee, and Number could be a shift. But let’s keep this abstract for now.

Say we want to assign each Letter to a Number, but it can only be assigned once. Now we could do this simply with a zip() operation, but we are going to make this problem subject to additonal constraints and rules later. Therefore, let’s do this mathematically and we will build on that foundation.

Envision a grid much like that board game Battleship where one axis is each Letter and the other is each Number. The intersection of each Letter or Number will be a 1 or 0 (true or false), indicating whether that Letter was assigned to that Number. Here is one possible solution done visually:

ONE TWO THREE FOUR FIVE
A 0 1 0 0 0
B 0 0 0 1 0
C 0 0 1 0 0
D 1 0 0 0 0
E 0 0 0 0 1

Bear with me. Let’s represent each “intersection” as a Slot class, which holds that given Letter and Number. We will also create an ojAlgo Variable that is binary, meaning it is constrained to be 1 or 0. We will also create a companion object to build and hold all the instances shown in this grid.

data class Slot(val letter: Letter, val number: Number) {
 
    val occupied = variable().binary()
 
    companion object {
        val all = Letter.values().asSequence().flatMap { letter ->
            Number.values().asSequence().map { number -> Slot(letter, number) }
        }.toList()
    }
    override fun toString() = "$letter$number: ${occupied?.value?.toInt()}"
}

How do you mathematically enforce a Letter only getting assigned to a Number once, and vice versa? Easy! To see if a letter has only been assigned once, just sum all of its “slot” values and enforce that it must not be greater than 1. Visually, all the values in a given column or a given row must sum to no more than 1. Mathematically, this would be expressed like this:

We could also use = instead of <=, but if we ever had more numbers than letters perhaps we would be okay with no assignment happening at all. We just do not want it assigned more than once.

Do the same for each Number. For instance, ONE's slots must sum to no more than 1.

To execute this simple logic, modify the Letter and Number classes to retrieve their slots. Then we will append an addConstraints() function to each one, where we will call addExpression(), set the upper() to 1, and then loop through each slot to “add” its variable. Since we have no need to multiply each variable with a coefficient, just set it to 1 as multiplying 1 will have no effect.

enum class Letter {
    A,
    B,
    C,
    D,
    E;
 
    val slots by lazy {
        Slot.all.filter { it.letter == this }.sortedBy { it.number.value }
    }
 
    fun addConstraints() {
 
        // constrain each letter to only be assigned once
        addExpression().upper(1).apply {
            slots.forEach {
                set(it.occupied, 1) // multiplier of 1
            }
        }
    }
}
 
enum class Number(val value: Int)  {
    ONE(1),
    TWO(2),
    THREE(3),
    FOUR(4),
    FIVE(5);
 
    val slots by lazy {
        Slot.all.filter { it.number == this }.sortedBy { it.letter }
    }
 
    fun addConstraints() {
 
        // constrain each NUMBER to only be assigned one slot
        addExpression().upper(1).apply {
            slots.forEach {
                set(it.occupied, 1) //multiplier of 1
            }
        }
    }
 
    override fun toString() = value.toString()
}

If that set() function is confusing you, we are basically just multiplying each binary by 1 which has no effect.

Now let’s run this. We will loop through each Letter and Number to call its addConstraints() function. We are not trying to minimize() or maximize() anything (such as profit or cost) in this problem. We just want something that is operable and works with all constraints we defined. Then we will print a pretty grid of all our results.

fun main(args: Array<String>) {
 
    Letter.values().forEach { it.addConstraints() }
    Number.values().forEach { it.addConstraints() }
 
    model.minimise().run(::println)
 
    Letter.values().joinToString(prefix = "   ", separator = "   ").run(::println)
 
    Number.values().forEach { n ->
        Letter.values().asSequence().map { l -> l.slots.first { it.number == n }.occupied.value.toInt() }
                .joinToString(prefix = "$n  ", separator = "   ").run { println(this) }
    }
}

Here is the output I got:

OPTIMAL 0.0 @ [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
   A   B   C   D   E
1  0   0   0   1   0
2  0   0   1   0   0
3  0   1   0   0   0
4  1   0   0   0   0
5  0   0   0   0   1

As you can see, we got an OPTIMAL result meaning it found a solution (as opposed to an INFEASIBLE if no solution is possible with our constraints). As you can see in our improvised console grid, each letter was cleanly assigned (via a 1) to a Number.

Now let’s jump this up a notch.

Restricting Items to Specific Regions

Say two additional rules need to be added to our model:

  • C must be assigned to a Number greater than or equal to THREE
  • D must be assigned to a Number less than or equal to TWO

These fixed constraints are usually not hard to implement. It is when we start making variables relative to other variables that things start to get gnarly (covered later in this article). But for these two rules, all we have to do is tease out the slots in those desired ranges, and say they must sum to 1.

We can add these two constraints to our model, and at least one of those slots for C and D respectively must be 1 for those equations to yield true. Let’s just do this in our main() function. We are going to use level() instead of upper() this to express a strict = rather than a <=.

fun main(args: Array<String>) {
 
    Letter.values().forEach { it.addConstraints() }
    Number.values().forEach { it.addConstraints() }
 
    // C must be greater than or equal to THREE
    addExpression().level(1).apply {
        Letter.C.slots.asSequence().filter { it.number.value >= 3 }.forEach {
            set(it.occupied, 1)
        }
    }
 
    // D must be less than or equal to TWO
    addExpression().level(1).apply {
        Letter.D.slots.asSequence().filter { it.number.value <= 2 }.forEach {
            set(it.occupied, 1)
        }
    }
 
 
    model.minimise().run(::println)
 
    Letter.values().joinToString(prefix = "   ", separator = "   ").run(::println)
 
    Number.values().forEach { n ->
        Letter.values().asSequence().map { l -> l.slots.first { it.number == n }.occupied.value.toInt() }
                .joinToString(prefix = "$n  ", separator = "   ").run { println(this) }
    }
}

When we run this, sure enough we get an output that will conform to those two new rules:

OPTIMAL 0.0 @ [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
   A   B   C   D   E
1  1   0   0   0   0
2  0   0   0   1   0
3  0   0   1   0   0
4  0   1   0   0   0
5  0   0   0   0   1

Alright, let’s escalate this into something more realistic.

Using Multiple Consecutive Blocks

When you start dealing with relative constraints, this is where things become less intuitive. Linear and integer programming limits you to only express functions in a simple addition/subtraction structure such as Ax + By + Cz <= 40 (A, B, and C are constants multiplied to variables x, y, and z). You cannot even do simple operations like multiply variables to other variables or use if expressions in the model functions. Therefore, you really have to challenge yourself to think creatively with pencil and paper, and not be afraid to experiment with mathematical equations and proofs to solve your problem before you even write code.

For instance, let’s say we want to add a concept to our model where a Letter can consume multiple contiguous slots. This is realistic as classrooms, staff, and other scheduling problems may break shifts up into say, 15 minute blocks, and a classroom can be scheduled for 90 minutes (or 6 blocks). In this case, let’s add the following rules:

  • Letter B requires 2 contiguous slots
  • Letter C requires 3 contiguous slots

The first thing we will need to do is add more Number items, as we do not have enough for these new requirements. Let’s add all the way up to NINE, even if it is actually more than we need (just because).

enum class Number(val value: Int)  {
    ONE(1),
    TWO(2),
    THREE(3),
    FOUR(4),
    FIVE(5),
    SIX(6),
    SEVEN(7),
    EIGHT(8),
    NINE(9);
 
    val slots by lazy {
        Slot.all.filter { it.number == this }.sortedBy { it.letter }
    }
 
    fun addConstraints() {
 
        // constrain each NUMBER to only be assigned one slot
        addExpression().upper(1).apply {
            slots.forEach {
                set(it.occupied, 1)
            }
        }
    }
 
    override fun toString() = value.toString()
}

So how in the world do we express this logic of multiple contiguous blocks as a linear function? It is not intuitive so the best place to start is to try and make it intuitive. Let’s plot with pencil and paper some valid scenarios for what we want. Here is an example of a desired scenario:

ONE TWO THREE FOUR FIVE SIX SEVEN EIGHT NINE
A 0 1 0 0 0 0 0 0 0
B 0 0 0 1 1 0 0 0 0
C 0 0 0 0 0 1 1 1 0
D 1 0 0 0 0 0 0 0 0
E 0 0 0 0 0 0 0 0 1

Let’s focus just on the row for C. Since C requires three contiguous slots, let’s break down a given case into rolling groups of three like this:

ONE TWO THREE FOUR FIVE SIX SEVEN EIGHT NINE
ALL 0 0 0 0 0 1 1 1 0
GROUP 1 0 0 0
GROUP 2 0 0 0
GROUP 3 0 0 0
GROUP 4 0 0 1
GROUP 5 0 1 1
GROUP 6 1 1 1
GROUP 7 1 1 0

Now we have rolling groups of three slots. GROUP 6 has the desired set of three contiguous 1’s. Let’s label each of these groups with a RESULT column, where 1 indicates a success and 0 indicates a fail like this:

RESULT ONE TWO THREE FOUR FIVE SIX SEVEN EIGHT NINE
ALL 0 0 0 0 0 1 1 1 0
GROUP 1 0 0 0
GROUP 2 0 0 0 0
GROUP 3 0 0 0 0
GROUP 4 0 0 0 1
GROUP 5 0 0 1 1
GROUP 6 1 1 1 1
GROUP 7 0 1 1 0

So what we really want is a series of constraints expressed as these equations:

Well that escalated quickly. Now before you throw in the towel, really study what is happening above and give it a chance to sink in. Pull out a pencil and paper if you need to and proof a few cases. Again, we broke C’s slots into contiguous, rolling groups of three. One group needs to result in a 1 while the rest must be a 0. The series of functions above express this. We create a new binary variable for each group that can be 1 or 0. We then add all those binaries together which should sum to 1, meaning only one of them must be 1 while the rest are 0.

The M is actually going to be 3 in this case, because we need three slots. We multiply it to the group’s binary and subtract it as a subtle way of executing a constraint. That expression must be greater than or equal to 0, and that M may or may not be voided if the group’s binary is 1 or 0. If we have 1’s scattered across multiple groups, the entire system of equations will invalidate that case. If three contiguous 1’s exist in a group, that will satisfy the system of equations.

Here is how we express this in our Kotlin model:

enum class Letter(val slotsNeeded: Int = 1) {
    A,
    B(slotsNeeded = 2),
    C,
    D(slotsNeeded = 3),
    E;
 
    val slots by lazy {
        Slot.all.filter { it.letter == this }.sortedBy { it.number.value }
    }
 
    fun addConstraints() {
 
        // constrain each LETTER to number of slots needed
        addExpression().level(slotsNeeded).apply {
            slots.forEach {
                set(it.occupied, 1)
            }
        }
 
        //only inflict model with contiguous logic if needed
        if (slotsNeeded > 1) {
 
            val allGroupSlots = addExpression().level(1)
 
            slots.rollingBatches(slotsNeeded).forEach { group ->
 
                val slotForGroup = variable().binary()
 
                allGroupSlots.set(slotForGroup, 1)
 
                addExpression().lower(0).apply {
                    group.forEach {
                        set(it.occupied,1)
                    }
                    set(slotForGroup, -1 * slotsNeeded)
                }
            }
        }
    }
}

When you run this, sure enough you get contiguous sets of slots occupied for B and C.

OPTIMAL 0.0 @ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
   A   B   C   D   E
1  0   0   0   0   0
2  0   0   0   1   0
3  0   0   0   1   0
4  0   0   0   1   0
5  0   0   1   0   0
6  0   1   0   0   0
7  0   1   0   0   0
8  1   0   0   0   0
9  0   0   0   0   1

Here is all the code in its entirety. You can also clone the code from GitHub.

import org.ojalgo.optimisation.ExpressionsBasedModel
import org.ojalgo.optimisation.Variable
import java.util.concurrent.atomic.AtomicInteger
 
// declare ojAlgo Model
val model = ExpressionsBasedModel()
 
// custom DSL for Example3.getModel expression inputs, eliminate naming and adding
val funcId = AtomicInteger(0)
val variableId = AtomicInteger(0)
fun variable() = Variable(variableId.incrementAndGet().toString().let { "Variable$it" }).apply(model::addVariable)
fun addExpression() = funcId.incrementAndGet().let { "Func$it"}.let { model.addExpression(it) }
 
 
fun main(args: Array<String>) {
 
    Letter.values().forEach { it.addConstraints() }
    Number.values().forEach { it.addConstraints() }
 
    // C must be greater than or equal to THREE
    addExpression().level(1).apply {
        Letter.C.slots.asSequence().filter { it.number.value >= 3 }.forEach {
            set(it.occupied, 1)
        }
    }
 
    // D must be less than or equal to TWO
    addExpression().level(1).apply {
        Letter.D.slots.asSequence().filter { it.number.value <= 2 }.forEach {
            set(it.occupied, 1)
        }
    }
 
 
    model.minimise().run(::println)
 
    Letter.values().joinToString(prefix = "   ", separator = "   ").run(::println)
 
    Number.values().forEach { n ->
        Letter.values().asSequence().map { l -> l.slots.first { it.number == n }.occupied.value.toInt() }
                .joinToString(prefix = "$n  ", separator = "   ").run { println(this) }
    }
}
 
enum class Letter(val slotsNeeded: Int = 1) {
    A,
    B(slotsNeeded = 2),
    C,
    D(slotsNeeded = 3),
    E;
 
    val slots by lazy {
        Slot.all.filter { it.letter == this }.sortedBy { it.number.value }
    }
 
    fun addConstraints() {
 
        // constrain each LETTER to number of slots needed
        addExpression().level(slotsNeeded).apply {
            slots.forEach {
                set(it.occupied, 1)
            }
        }
 
        // constrain slots to be consecutive if more than one are needed
        // this is tricky because slots need to be "batched" in rolling groups of needed slot size
        // e.g. for a slot size three, we need (x1,x2,x3), (x2,x3,x4), (x3,x4,x5) and so on. A binary is attached to each group
        // and another binary needs to be shared across all the batches
        // x1 + x2 + .. xi - Sb >= 0
        // Where xi is slot binaries, S is number of slots needed, and b is shared binary across all groups
        if (slotsNeeded > 1) {
 
            val allGroupSlots = addExpression().level(1)
 
            slots.rollingBatches(slotsNeeded).forEach { group ->
 
                val slotForGroup = variable().binary()
 
                allGroupSlots.set(slotForGroup, 1)
 
                addExpression().lower(0).apply {
                    group.forEach {
                        set(it.occupied,1)
                    }
                    set(slotForGroup, -1 * slotsNeeded)
                }
            }
        }
    }
}
 
enum class Number(val value: Int)  {
    ONE(1),
    TWO(2),
    THREE(3),
    FOUR(4),
    FIVE(5),
    SIX(6),
    SEVEN(7),
    EIGHT(8),
    NINE(9);
 
    val slots by lazy {
        Slot.all.filter { it.number == this }.sortedBy { it.letter }
    }
 
    fun addConstraints() {
 
        // constrain each NUMBER to only be assigned one slot
        addExpression().upper(1).apply {
            slots.forEach {
                set(it.occupied, 1)
            }
        }
    }
 
    override fun toString() = value.toString()
}
 
data class Slot(val letter: Letter, val number: Number) {
 
    val occupied = variable().binary()
 
    companion object {
        val all = Letter.values().asSequence().flatMap { letter ->
            Number.values().asSequence().map { number -> Slot(letter, number) }
        }.toList()
    }
    override fun toString() = "$letter$number: ${occupied?.value?.toInt()}"
}
 
fun <T> List<T>.rollingBatches(batchSize: Int) = (0..size).asSequence().map { i ->
    subList(i, (i + batchSize).let { if (it > size) size else it })
}.filter { it.size == batchSize }

Conclusions

Linear programming is fairly easy to model in most cases, but as you can see here integer/binary programming can be hard. I am still learning binary patterns to this day, and find it helpful to spend a lot of time with pencil and paper experimenting with equations. I find it is a hit-or-miss process with lots of proof work, and I hope to get proficient in tackling even larger real-world problems.

Thankfully Kotlin made it easy, organized, and refactorable to execute, and allowed us to spend more time thinking about the math. As you probably have noticed, the Kotlin stdlib came incredibly handy in preparing data for the model.

If you are dealing with thousands of variables, you may encounter performance issues even with ojAlgo. Thankfully most solvers, including ojAlgo, support plugging in commercial optimizer implementations like CPLEX. When your solver is acting exhausted or slow, and you are sure there is no modeling mistake, it may be worthwhile getting a commercial-grade solver.

In the coming few months, I will do a few more posts on using Kotlin for linear/integer programming and optimization. Until then, be sure to check out this book as well as this one. Also be sure to use Math Stack Exchange for specific modeling questions, and StackOverflow for ojAlgo questions.

No comments:

Post a Comment