lpSolve with absolute values

    ntypes <- length(unique(dat[,"type"])) # number of types
    typemap <- setNames(seq.int(ntypes), unique(dat[,"type"])) # map typename to 1,...,ntypes

    solve_one <- function(subdat, capdat) {

      # create object
      lprec <- make.lp(0, ncol=3*ntypes) # for each type, two decision variables, one for even-ness

      # By convention, we assume that the first ntypes variables are physics for type 1, ..., ntypes
      # and the second ntypes variables are maths

      # add objective and type
      set.objfn(lprec, obj=c(rep(100, ntypes), rep(65, ntypes), rep(1, ntypes)))

      # add capacity constraints
      idx <- which(capdat[,"mon"]==subdat[1,"mon"] & capdat[,"region"]==subdat[1,"region"]) # lookup the right cap
      add.constraint(lprec, rep(1, ntypes), type="<=", rhs=capdat[idx,"regutil"], indices=seq.int(ntypes))
      add.constraint(lprec, rep(1, ntypes), type="<=", rhs=capdat[idx,"nonregutil"], indices=seq.int(ntypes+1, 2*ntypes))

      # add allsub equality constraints and minimum constraints
      for (typ in subdat[,"type"]) {
        add.constraint(lprec, c(1,1), type="=", rhs=subdat[typemap[typ], "effort_calculated"], indices=c(typemap[typ], ntypes+typemap[typ]))
        add.constraint(lprec, 1, type=">=", rhs=subdat[typemap[typ],"regmin"], indices=typemap[typ])
        add.constraint(lprec, 1, type=">=", rhs=subdat[typemap[typ],"nonregmin"], indices=ntypes+typemap[typ])
      }

      # set last decision vars as absolute difference between nonreg and reg for each type
      # see here for an explanation of the trick: http://lpsolve.sourceforge.net/5.1/absolute.htm
      for (i in seq.int(ntypes)) {
        add.constraint(lprec, c(1,-1,-1), type="<=", rhs=0, indices=c(i, ntypes+i, 2*ntypes+i))
        add.constraint(lprec, c(-1,1,-1), type="<=", rhs=0, indices=c(i, ntypes+i, 2*ntypes+i))
      }

      # solution data.frame
      ans <- subdat[, c("mon", "region", "type")]

      # solve
      if(solve(lprec)==0) {
        sol <- get.variables(lprec)
        for (i in seq.int(nrow(subdat))) {
          ans[i, "regval"] <- sol[typemap[subdat[i,"type"]]]
          ans[i, "nonregval"] <- sol[typemap[subdat[i,"type"]]+ntypes]
        }
      } else ans[,c("regval", "nonregval")] <- NA # no solution found

      return(ans)
    }

    sp <- split(dat, list(dat[,"mon"], dat[,"region"]))
    # browser()x
    results <- lapply(sp, solve_one, capdat=capdat)
    results <- do.call(rbind, results)
    rownames(results) <- NULL
    results