This neat approach showed up recently as an answer to a FiveThirtyEight puzzle and of course I couldn’t help but throw it at dplyr as soon as I could. Turns out that’s not a terrible idea. The question posed is
optimise
200a + 100b + 50c + 25d
under the constraints
400a + 400b + 150c + 50d ≤ 1000, b ≤ a, a ≤ 1, c ≤ 8, d ≤ 4,
and (a,b,c,d) all non-negative integers.
Leaving aside any interpretations of wording of the original question (let’s just start with trying to solve this system of inequalities) the solution provided used 4 nested loops, which can definitely be avoided.
My approach was to create all possible combinations of the 4 variables (within the given constraints), filter out the ones that don’t meet the constraint criteria, then sort by the evaluating expression to find which one does best.
I’m not suggesting that this is by any means always the best approach, but when the phase-space of possible solutions is so low (especially combinations of small integers) then this is pretty tidy (technically a single dplyr chain).
Alternatively, one could set this up as an equation and use a linear solver. In that case, we want to optimise \(\max(\|A x\|)\) subject to the constraints \(G x \ge h\) where \(A\) and \(x\) represent the coefficients and variables to be optimised, \(h\) the constraint vector, and \(G\) a matrix of coefficients for the constraints. For the system we’re looking at, that matrix inequality looks like this
\[ \left[\begin{array}{cccc}400 & 400 & 150 & 50 \\1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{array}\right] \left[\begin{array}{c}a \\ b \\ c \\ d\end{array}\right] \le \left[\begin{array}{c}1000 \\ 1 \\ 1 \\ 8 \\ 4\end{array}\right]\ . \]
Of course, the constraint that \(b \le a\) needs to be checked after the fact.
Programming this is fairly straightforward, even with the constraint that these are integer solutions. limSolve::linp is made for exactly these types of problems.
which results in the same answer as our manual brute-force search.
One last thing to try is to plot the solution space and see how it looks. Sounds like a good opportunity to try out plotly.
Since this is technically a 5D plot (4 variables and a value) it’s a little difficult to visualise. I’ve reduced the dimensionality by treating each unique combination of \(a \le 1\) and \(b \le a\) (i.e. \(00,~01,~10,~11\)) as a group and using colour to distinguish those. The plot below should show up as a 3D object, so click, drag, and scroll it to have a closer look. Clicking on a group will remove/add it so you can get a clearer view, and hovering over a point should bring up the values of the axes and evaluation.
Going back to the expression that’s being optimised it’s pretty clear why it’s broken down into 4 planes when grouped this way (substitute different values of \(a\) and \(b\) to see).
Do you have another way to solve this? Drop a line or a link in the comments.