*Towards Data Science*.

## Friday, January 19, 2018

## Saturday, January 13, 2018

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