Kotlin and Linear Programming Part II - Linear Optimization

In my previous article in this series, I introduced using linear programming with Kotlin. Again, linear programming broadly encompasses linear, integer, and binary programming all of which can be helpful for optimization problems. It was probably unconventional I started with binary programming, but it arguably makes linear programming useful beyond simple mixing problems. For instance, binary programming enables tackling complex scheduling problems which I will revisit again in my next post. But for now, let’s get pure linear programming out of the way. Then I will get to the exciting stuff.

Linear Programming and Optimization

Linear programming, in its purest form, optimizes continuous variables. Continuous means a variable can optimize to a decimal precision rather than a whole number. Unlike integer programming where a variable evaluates to a whole integer like 2, linear programming can optimize it to 2.4, 2.4012 or even 2.04000125. Surprisingly, this is much easier for a machine to solve than integer or binary programming. The machine can gracefully solve variables that truly are linear, rather than rely on brute-force approximations that integer programming requires. When you mix linear and integer programming together, it is called mixed programming.

One thing we really did not touch on in the last article is optimization. Our abstract binary problem was not given a goal other than find a solution that was “feasible” and satisfied the constraints. But in optimization, many solutions can exist. For example, an optimizer can produce many possible manufacturing plans (produce this many items of Product X, and this many Product Y), but it is likely we want the one that is the least costly, or most profitable. Therefore we can provide a “goal” function that is used to choose the best solution. This function will mathematically manipulate the variables in a way that expresses cost, revenue, machine/staff utilization, etc and the algorithm will minimize or maximize that function. You can even have multiple goals put into one function, with each one weighted based on their importance.

The Driver Problem

There are a lot of linear problem examples out there that deal with blending, from manufacturing cars to making sausage. But I am interested in scheduling at the moment so let’s put a linear twist on that, even if it conventionally is an integer problem.

At KotlinConf, I presented a driver scheduling example which garnered interest. This is a linear problem with continuous and binary variables, therefore making it a mixed problem. Driver shift start and end variables are allowed to move freely and continuously throughout the day, and the shift must be 4-6 hours in length. This continuous definition of schedule shifts is not likely how schedules are built in real life (they are usually discrete and integer-based), but it is an interesting application nonetheless. Let’s take a look at it:

You have three drivers who charge the following rates:

Driver 1: $10 / hr
Driver 2: $12 / hr
Driver 3: $15 / hr

From 6:00 to 22:00, schedule one driver at a time to provide coverage, and minimize cost.

Each driver must work 4-6 hours a day. Driver 2 cannot work after 11:00.

We could expand the scope of this problem significantly, scheduling an entire week and putting additional constraints on that scope (e.g. a worker cannot be scheduled more than 30 hours/week). But let’s keep the scope of this problem limited to one day for now, and save multi-day scheduling for the next article.

Structuring the Problem

Let’s start by saving the operating day and allowable shift size as Kotlin ranges. If we represent the day as 24 hours, the operating day will be 6:00 to 22:00, or 6..22 as a Kotlin IntRange. The allowableShiftSize can only be 4 to 6 hours, which would be an IntRange of 4..6.

val operatingDay = 6..22
val allowableShiftSize = 4..6

Declare an instance of an ExpressionBasedModel. This is the entity that we will input constraint functions into. Let’s also define an improvised Kotlin DSL to streamline ojAlgo. These variable() and addExpression() functions will automatically assign a name to each variable and function, which personally I find tedious to do. Let’s also calculate the operatingDayLength and save it to a constant for convenience.

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) }


// constants
val operatingDayLength = operatingDay.endInclusive - operatingDay.start

Now let’s get to the core of the problem.

data class Driver(val driverNumber: Int,
                  val rate: Double,
                  val availability: IntRange? = null) {

    val shiftStart = variable().lower(6).upper(22)
    val shiftEnd = variable().lower(6).upper(22)

    fun addToModel() {
        // expression inputs go here
    }
}

As shown above, create a Driver class with two ojAlgo variables, which will contain the shiftStart and shiftEnd values after the model is optimized. Note these two variables are bound to a lower() of 6 and an upper() of 22 which is our operating day. Also declare a driverNumber which will simply act as an ID, the rate which will be type Double, and a nullable IntRange called availability for drivers that can work only certain hours. Finally, put an empty addToModel() function which is where we will write code that “puts” the driver mathematically into the ojAlgo model. We will get to this in a moment.

We can then hold a list of Drivers. As a placeholder, we will have a main() function that calls a buildModel() function, then minimize() the cost, and print the results.

// declare drivers
val drivers = listOf(
        Driver(driverNumber = 1, rate = 10.0),
        Driver(driverNumber = 2, rate = 12.0, availability = 6..11),
        Driver(driverNumber = 3, rate = 14.0)
)


// parameters
val operatingDay = 6..22
val allowableShiftSize = 4..6

fun main(args: Array<String>) {

    buildModel()

    model.minimise().run(::println)

    // see variables for each driver
    drivers.forEach {
        println("Driver ${it.driverNumber}: ${it.shiftStart.value.toInt()}-${it.shiftEnd.value.toInt()}")
    }
}

fun buildModel() {

}

Implementing the Math

Next we need to populate the buildModel() function with our mathematical constraints. Some of those constraints will be specific to each Driver, in which case those will be in the Driver’s addToModel() function. Other constraints will apply to all the drivers.

Day Coverage

We have already constrained the shiftStart and shiftEnd variables to the operating day, so we are good there (using the lower() and upper() functions). We still need a constraint that ensures our entire 16-hour operating day has been covered. This is pretty easy to express mathematically. The sum of differences between each shiftEnd and shiftStart for each driver must sum up to the operating day’s length:

Therefore in our buildModel() function we can call addExpression(), level() it to the operatingDayLength (which should be 16), and use a handy Kotlin apply() closure to loop through each driver. We add each driver’s shiftStart and shiftEnd variables to the expression with the appropriate 1 or -1 multiplier to add or subtract them.

fun buildModel() {

    //ensure coverage of entire day
    model.addExpression()
            .level(operatingDayLength)
            .apply {
                drivers.forEach {
                    set(it.shiftEnd, 1)
                    set(it.shiftStart, -1)
                }
            }
}

Cost Objective

Before we create expressions in context with each individal Driver, let’s define our objective(). We want to minimize cost, but for now we just need an expression that calculates the cost which we call the minimise() function against. This expression will be similar to the day coverage one, but we multiply each difference by the driver’s rate. This will yield the total cost.

Then we express this in the buildModel() function in a similar way as the day coverage, but when we call addExpression() and give it a weight(1) of 1. This will make the expression a goal driver rather than a constraint. Note we can actually provide multiple objectives and give each one a weight based on their importance, but we will stick with 1 for now.

fun buildModel() {

    //ensure coverage of entire day
    model.addExpression()
            .level(operatingDayLength)
            .apply {
                drivers.forEach {
                    set(it.shiftEnd, 1)
                    set(it.shiftStart, -1)
                }
            }

    // set objective
    model.expression().apply {
        weight(1)
        drivers.forEach {
            set(it.shiftEnd, it.rate)
            set(it.shiftStart, -1 * it.rate)
        }
    }

    // driver-specific expressions
    drivers.forEach { it.addToModel() }
}

Let’s also add a call to each driver’s addToModel() function which we will implment next.

Shift Length

A driver must work a minimum of 4 hours but a maximum of 6 hours. This one is pretty simple. We just express the difference of the shiftEnd and shiftStart as an inequality.

In the Driver’s addToModel() function, we can call addExpression(), constrain its lower() and upper() bounds to the allowable shift size range. Then we set the shiftStart and shiftEnd variables, multiplying the shiftStart by -1 so it is subtracted rather than added.

data class Driver(val driverNumber: Int,
                  val rate: Double,
                  val availability: IntRange? = null) {

    val shiftStart = variable().lower(6).upper(22)
    val shiftEnd = variable().lower(6).upper(22)

    fun addToModel() {

        //constrain shift length
        model.addExpression()
                .lower(allowableShiftSize.start)
                .upper(allowableShiftSize.endInclusive)
                .set(shiftEnd, 1)
                .set(shiftStart, -1)

    }
}

Shift Restrictions

Driver 2 is the only driver that restricted their shift to something less than the operating day. He wants to only work within the 6 to 11 range. This constraint is expressed as two very simple inequalities:

If a Driver indeed has provided a specific availability that is not null, we can take that range and use it to build those constraints as shown below:

data class Driver ... {

    fun addToModel() {

        //constrain shift length
        model.addExpression()
                .lower(allowableShiftSize.start)
                .upper(allowableShiftSize.endInclusive)
                .set(shiftEnd, 1)
                .set(shiftStart, -1)

        //add specific driver availability
        availability?.let {
            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftStart, 1)

            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftEnd, 1)
        }
    }
}

Preventing Shift Overlap

Okay, everything else up to this point was fairly easy, right? Now we come to the hard part and introduce some binary variables to prevent shift overlap. If you run with just these constraints we put in so far, we are going to have drivers overlapping their shifts! We need to prevent this and unfortunately, we have to think very abstractly to do it.

To prevent overlaps, each driver needs to have a mathematical relationship to the other drivers. Let’s focus on Driver 1 and Driver 2. The shiftStart of Driver 1 may be after the shiftEnd of Driver 2. But the shiftStart of Driver 2 could also be after the shiftEnd of Driver 1. Both cases are visually shown below:

  S2________E2          S1__________E1

  S1________E1          S2__________E2

Both are valid cases, and we need either to be true as a constraint to prevent overlap. But how do we support both cases in our mathematical model, which expects all constraints to be satisfied in order to yield a solution? You logically cannot declare both cases must be true, right? It seems like a paradox. We need to an “OR” rather than an “AND”.

So how do we mathematically express that “OR”? There is a way! For a given relationship between a driver i and another driver j, you can share a boolean variable between them that must be 1 or 0. Here is a modification to our constraints that effectively achieves the “OR”:

Well, that jumped up a notch. The M is going to be the length of the planning window, which is 16 hours. We just need M to be a very large number, so we could make it a 1000 but the window size should be fine.

The binary variable is represented by a delta, and it can only be optimized to 1 or 0. Notice the subtle math this achieves and allows both cases of Si coming before Ej, or Sj coming before Ei, but demands at least one of them must be true. M is a large enough number (the planning window) to throw the inequality far out enough that it allows the “false” case to now be true. It effectively nullifies. Thus, this achieves our much-needed “OR” operation for our constraint.

If your head is spinning pull out a pencil and paper, throw a few values at the expression, and you’ll see what I mean.

Here is how we can implement this for our Driver instance. It must establish this constraint with each other Driver. Note that ojAlgo needs all variables on one side of the equation and distributed, so we may have to do a little bit of basic algebra to move our variables around like this:

data class Driver ... {

    fun addToModel() {

        //constrain shift length
        model.addExpression()
                .lower(allowableShiftSize.start)
                .upper(allowableShiftSize.endInclusive)
                .set(shiftEnd, 1)
                .set(shiftStart, -1)

        //add specific driver availability
        availability?.let {
            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftStart, 1)

            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftEnd, 1)
        }

        //prevent shift overlap
       drivers.asSequence()
               .filter { it != this }
               .forEach { otherDriver ->

                   val occupied = variable().binary()

                   model.addExpression()
                           .upper(0)
                           .set(otherDriver.shiftEnd, 1)
                           .set(occupied, operatingDayLength * - 1)
                           .set(shiftStart, -1)

                   model.addExpression()
                           .upper(operatingDayLength)
                           .set(shiftEnd, 1)
                           .set(occupied, operatingDayLength)
                           .set(otherDriver.shiftStart, -1)
               }
    }
}

And there you have it. If you find this is dizzying to grasp, spend some time with pencil and paper, visualizing how this equation works and throw different scenarios and variable values at it. I have found visualizing the problem is the most effective way to make it intuitive.

Running the Program

Finally, go back to your main() function and call the minimise() function, and print the results of the variables like this:

fun main(args: Array<String>) {

    buildModel()

    model.minimise().run(::println)

    // see variables for each driver
    drivers.forEach {
        println("Driver ${it.driverNumber}: ${it.shiftStart.value.toInt()}-${it.shiftEnd.value.toInt()}")
    }
}

You should get the following output, with these shifts that will minimize the cost according to your constraints.

OPTIMAL 190.0 @ [11.0, 17.0, 6.0, 11.0, 17.0, 22.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0]
Driver 1: 11-17
Driver 2: 6-11
Driver 3: 17-22

Hopefully you found this somewhat interesting and exciting, and this topic of linear/integer programming is giving you new ways to look at optimization problems. In the next article, I will cover an ambitious scheduling example that is more real-life, and schedules classrooms and classes across an entire week.

Here is the entire code example. You can also view it on GitHub here.

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 model 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 ExpressionsBasedModel.addExpression() = funcId.incrementAndGet().let { "Func$it"}.let { addExpression(it) }

// parameters
val operatingDay = 6..22
val allowableShiftSize = 4..6

// constants
val operatingDayLength = operatingDay.endInclusive - operatingDay.start


// declare drivers
val drivers = listOf(
        Driver(driverNumber = 1, rate = 10.0),
        Driver(driverNumber = 2, rate = 12.0, availability = 6..11),
        Driver(driverNumber = 3, rate = 14.0)
)


fun main(args: Array<String>) {

    buildModel()

    model.minimise().run(::println)

    // see variables for each driver
    drivers.forEach {
        println("Driver ${it.driverNumber}: ${it.shiftStart.value.toInt()}-${it.shiftEnd.value.toInt()}")
    }
}

fun buildModel() {

    //ensure coverage of entire day
    model.addExpression()
            .level(operatingDayLength)
            .apply {
                drivers.forEach {
                    set(it.shiftEnd, 1)
                    set(it.shiftStart, -1)
                }
            }

    // set objective
    model.objective().apply {
        drivers.forEach {
            set(it.shiftEnd, it.rate)
            set(it.shiftStart, -1 * it.rate)
        }
    }

    // driver-specific expressions
    drivers.forEach { it.addToModel() }
}

// Driver class will put itself into the Model when addToModel() is called
data class Driver(val driverNumber: Int,
                  val rate: Double,
                  val availability: IntRange? = null) {

    val shiftStart = variable().lower(6).upper(22)
    val shiftEnd = variable().lower(6).upper(22)

    fun addToModel() {

        //constrain shift length
        model.addExpression()
                .lower(allowableShiftSize.start)
                .upper(allowableShiftSize.endInclusive)
                .set(shiftEnd, 1)
                .set(shiftStart, -1)

        //add specific driver availability
        availability?.let {
            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftStart, 1)

            model.addExpression()
                    .lower(it.start)
                    .upper(it.endInclusive)
                    .set(shiftEnd, 1)
        }

        //prevent shift overlap
        drivers.asSequence()
                .filter { it != this }
                .forEach { otherDriver ->

                    val occupied = variable().binary()

                    model.addExpression()
                            .upper(0)
                            .set(otherDriver.shiftEnd, 1)
                            .set(occupied, operatingDayLength * - 1)
                            .set(shiftStart, -1)

                    model.addExpression()
                            .upper(operatingDayLength)
                            .set(shiftEnd, 1)
                            .set(occupied, operatingDayLength)
                            .set(otherDriver.shiftStart, -1)
                }
    }
}