http://jcarroll.com.au/2016/04/27/solving-inequality-the-math-kind/

inequality_1.R

## 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=400*a + 400*b + 150*c + 50*d,
         optim.fun=200*a + 100*b +  50*c + 25*d) %>%
  filter(b <= a, const.fun <= 1000) %>%
  arrange(-optim.fun) %>%
  head(1)
## a b c d
## 1 0 3 3
## 
## optim.fun = 425

inequality_2.R

## 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