# Modeling state in linear programs

# Modeling state in linear programs

There’s a neat way to model decision making as an (integer) linear program for managing a collection of resources. In my head it’s called the “state-dynamics” model. In other places I’ve seen it, such as this set of notes, it’s an unnamed “formulette,” and they describe it as a “periodic review of a continuous process.” That name also makes sense, because you can view this model as periodically sampling a continuous function, and optimizing over an (approximate) integral of the function.

I’ll explain the situation with bananas. Say you’re managing an inventory of bananas. You occasionally get orders from customers to buy your bananas, and you occasionally place orders for new bananas from wholesale suppliers. As you read this, you may substitute “bananas” for any good that you can bucket over time and according to specific features.

You can model this by introducing “state” variables
representing the quantity of bananas in stock,
say `Quantity_t`

,
where `t`

is a regularly sampled time across
some predetermined planning horizon,
say monthly samples for a year.
Then you add “dynamics” constraints that look like

Quantity_{t+1} = Quantity_t + Purchase_[t,t+1] - Order_1 - Order_2 - ...

Where `Purchase_[t,t+1]`

is an integer variable
representing the amount of bananas
the solver will choose
for you to purchase during the month,
and `Order_i`

is the quantity ordered
by a customer order
(which, say, occurs once per order instance).

Then you add constraints like

Quantity_0 = CURRENT_STOCK Quantity_t >= SAFETY_FLOOR for all t

The solver gets to “solve” for the `Quantity_i`

,
but any system that respects the dynamics constraints
is literally tracking the state of banana inventory each month.
The dynamics constraints ensure that
each subsequent state variable is consistent
with the previous state and the changes in each interval.
And the `CURRENT_STOCK`

constraint is the base case.
So the only “choices” for the state variables
are forced through the *real* decision variables,
which are the purchases.
Though, the solver doesn’t “know” that
when it solves the model.

The optimization objective can be something like, “minimize the sum of purchase costs plus the sum of inventory storage costs,” which can be written as a linear function of the purchase variables and state variables, respectively. Now you have yourself a tidy little inventory management model. Because we’re sampling over time, the total cost is more like an integral over time, and so you’d want to apply a quadrature rule to the samples to get area. A good one that is linear is the trapezoid rule.

This model seems quaint, but it has some hidden surprises that make it a very flexible and efficient formulation.

## Few integer variables

In solving an integer program, the performance is closely related to the number of integer variables the solver must solve for. This is because non-integer variables can be solved efficiently using the simplex method (and related methods), while integer variables must go through a branch-and-bound exponential search.

The state-dynamics model above is nice because
the only integer variables are the purchase variables.
That is, the state variables can all be marked as continuous,
and treated as continuous by the solver,
and it just so happens that if `CURRENT_STOCK`

is an integer,
and if each `Order_i`

is an integer,
then all the state variable values
are forced to be integer for free,
because integer values are the only feasible solutions.
The solver will find them without a fuss.

This is a commonly used trick, where the organization of the model ensures some variables must be integer in any feasible solution. If you can do this with ALL your variables, then you can use a much faster linear solver. In our case, the purchase variables still need to be integer.

## Easy to generalize

Imagine you have multiple brands of bananas now,
and you have to decide which bananas to use to fulfill
customer orders. In this case, the `Quantity_t`

state variables gain an index,
like `Quantity_{Dole, t}`

or `Quantity_{Chiquita, t}`

,
and the corresponding dynamics variables
and purchase variables get a new index.

Then the `Order_i`

become new (integer) decision variables
instead of just constants,
like `Order_{i, Dole}`

,
whose solution value counts
how many Dole bananas are used to fulfill order `i`

.

Similarly, your inventory might be spread
across many warehouses,
with time delays between purchasing some types of bananas
and having them arrive at a warehouse.
This can all be incorporated into the dynamics constraints,
with a little bit of extra work to track the relation
between a purchase in month `t_1`

and the later month `t_2`

that the purchased inventory arrives
and can be used to fulfill orders.

Another nice way the dynamics constraints generalize, is that they can incorporate any sort of fixed changes to the inventory stock. For example, you might have committed purchase orders for new bananas, and committed customer orders. These would be incorporated into the dynamics constraints as constants instead of variables.

## Easy to avoid infeasibility

A natural version of this model involves purchase limits, i.e., statements from your supplier that they can only deliver so much.

In such cases,
you might get into a situation
in which the banana demand
exceeds your supply at some point during the year.
Without any special work,
the solver will crash and say “this model is infeasible.”
In this case,
the constraints `Quantity_{Dole, t} >= SAFETY_FLOOR`

are impossible to satisfy,
because orders force it to dip below the floor,
but it is not possible to purchase enough bananas to recover,
due to the purchase limits or time delays.
Infeasibility grinds the system to a halt,
preventing the solver from even making partial
decisions about what to buy or fulfill.

A nice workaround is to use slack variables.
In this case, you’d add a new continuous variable
`z_{Dole, t} >= 0`

and modify the safety floor constraint to be
`Quantity_{Dole, t} + z_{Dole, t} >= SAFETY_FLOOR`

.
This allows the solver to always satisfy this constraint,
so to prevent the solver from using the `z`

unless absolutely necessary,
you add an objective term of
`+ C * z_{Dole, t}`

,
where `C`

is a sufficiently large constant,
say, 100x the largest non-slack coefficient in the objective.
This has the effect of making fake bananas
infinitely available for purchase,
but their “purchase cost” is too prohibitive
for the fake bananas to be optimal unless
the model is otherwise infeasible.

If every non dynamics constraint has this sort of slack variable incorporated, then the model can never be infeasible by design. On top of that, the slack variables act as an alerting mechanism. If a slack variable is nonzero in the solution, your system should report it to a human operator. This means that there is some defect that cannot be corrected without human intervention. Maybe they will cancel some orders, or rush-order more bananas to fix the problem. If the slack variables are set up well, each reported defect is localized (see next section), and the “blast radius” of the defect can be calculated.

## Bugs with alternative slack formulations

This slack variable business
is a lot more nuanced than it seems!
First, the choice of `C`

is tricky,
since if it’s too large it will cause numerical instability,
and if it’s too small the solver will use it to cheat
and invent bananas out of thin air.

Second, it’s not clear how critical it is
that I put the slack variable
in the safety floor constraint.
For example, you might think to put
the `z_{Dole,t}`

in the dynamics constraint.
It might seem especially good to do it that way
if the infeasibility appears to come from unclean data,
because it would “correct” the error and keep it corrected
in later months.

Doing it that way introduces subtle bugs, which I’ve personally made, and I’ve watched colleagues make but missed it in my review. The following situation explains. Say your initial banana inventory/order data is incorrect. In one warehouse, the starting inventory is 100 bananas, and the first month has 200 orders worth of bananas, which is infeasible. Maybe someone forgot to mark an order as canceled. Imagine further that all banana purchases take one month to arrive, so the solver can’t buy new bananas and is forced to use slack variables to stay feasible.

Now say you put the slack variable
in the dynamics constraint,
and there’s an incoming (committed) purchase order for 500 bananas
arriving the next month.
The constraints would look like (`Q`

for `Quantity`

)

0 <= Q_{Dole, Jan} = 100 0 <= Q_{Dole, Feb} = Q_{Dole, Jan} - 200 + z_{Dole, Jan/Feb} 0 <= Q_{Dole, Mar} = Q_{Dole, Feb} + 500 + z_{Dole, Feb/Mar} 0 <= Q_{Dole, Apr} = Q_{Dole, Mar} - Order_{1, Dole} + z_{Dole, Mar/Apr}

Assume the safety floor values are all zero.
The value of `z_{Dole, Jan/Feb}`

is forced to be 100,
but then the order of 500 the next month
makes the state variable in March be 500!
This means that the solver will think it has 500 Dole bananas
to use to fulfill Order 1, when 100 of those are faked.
Worse, it will be harder to trace where those fake bananas came from,
since it was a slack variable a few months earlier in the model.
The impact of correcting the defect
persisted beyond the period in which the deficit was present.

By contrast if the slack variable
is in the safety floor constraint,
then the value of `Q_{Dole, Feb}`

will be -100,
and when it comes time to fulfill Order 1,
we will only have 400 Dole bananas to pull from.
This is the right thing to do because
it’s the risk-averse approach.
We can’t know whether the 100 banana deficit
corresponds to a real order
(which might be pushed back a month)
or a data error,
so it’s better to debit it from the future 500 banana credit
and make the solver compensate
by using Chiquitas to partially fulfill Order 1,
rather than promise to dole out phantom Doles.

Pushing the slack variable to the safety floor constraint achieves this, and it further allows you to easily connect slack variable usage alerts to the orders they impact, since they have to be in the dynamics constraints in adjacent time periods.