Solving Inequality (the math kind) -

This neat approach (xianblog.wordpress.com) showed up recently as an answer to a FiveThirtyEight puzzle (fivethirtyeight.com) 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 (xianblog.wordpress.com), 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.

	## Optimising an expression subject to inequality constriaints
	## as seen on
	## http://fivethirtyeight.com/features/you-have-1-billion-to-win-a-space-race-go/
	## https://xianblog.wordpress.com/2016/04/21/an-integer-programming-riddle/
	## http://www.r-bloggers.com/an-integer-programming-riddle/
	## Blogged @ http://jcarroll.com.au/2016/04/27/solving-inequality-the-math-kind/

	## The expression to optimise:
	## 200a + 100b + 50c + 25d
	##
	## The constraints
	## 400a + 400b + 150c + 50d <= 1000
	## b <= a, a <= 1, c <= 8, d <= 4

	## load required packages
	library(dplyr)

	## optimise the expression by brute-forcing all integer
	## combinations, enforcing constraints, and sorting
	expand.grid(a=0:1, b=0:1, c=0:8, d=0:4) %>%
	mutate(const.fun=400a + 400b + 150c + 50d,
	optim.fun=200a + 100b + 50c + 25d) %>%
	filter(b <= a, const.fun <= 1000) %>%
	arrange(-optim.fun) %>%
	head(1)
	## a b c d
	## 1 0 3 3
	##
	## optim.fun = 425

view raw (gist.github.com) inequality_1.R (gist.github.com) hosted with ❤ by GitHub (github.com)

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 (www.inside-r.org) is made for exactly these types of problems.

	## Optimising an expression subject to inequality constriaints
	## as seen on
	## http://fivethirtyeight.com/features/you-have-1-billion-to-win-a-space-race-go/
	## https://xianblog.wordpress.com/2016/04/21/an-integer-programming-riddle/
	## http://www.r-bloggers.com/an-integer-programming-riddle/
	## Blogged @ http://jcarroll.com.au/2016/04/27/solving-inequality-the-math-kind/

	## The expression to optimise:
	## 200a + 100b + 50c + 25d
	##
	## The constraints
	## 400a + 400b + 150c + 50d <= 1000
	## b <= a, a <= 1, c <= 8, d <= 4

	## load required packages
	library(limSolve) ## limSolve::linp

	## alternatively, set the system up as a matrix
	## inequality constraints (LHS)
	## negative for <=
	G = -matrix(c(400,400,150,50,
	1,0,0,0,
	0,1,0,0,
	0,0,1,0,
	0,0,0,1), nrow=5, byrow=TRUE)

	## inequality constraints (RHS)
	## negative for <=
	H = -c(1000, 1, 1, 8, 4)

	## cost function (negative to maximise)
	Cost <- -c(200, 100, 50, 25)

	## optimise the system
	sol <- linp(Cost=Cost, G=G, H=H, int.vec=c(1,2,3,4))

	## the solution matches the one found in the dplyr chain
	sol$X
	## a b c d
	## 1 0 3 3
	-sol$solutionNorm
	## 425

view raw (gist.github.com) inequality_2.R (gist.github.com) hosted with ❤ by GitHub (github.com)

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 (plot.ly).

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.

rstats

See also