Distribute value to zones using population ratio and min/max criteria in R

This function carries out the algorithm you describe.

It first checks whether you have enough hats to make the minimum requirements. If not, it throws an error.

It then sees whether there are more than enough hats to go round, in which case it gives the maximum hats.

Otherwise, it distributes the minimum number of hats and subtracts this sum from remaining hats. It will then loop through, giving a hat to the remaining row with the largest gap between current hats and maximum hats multiplied by population size until there are no hats left to distribute.

distribute_hats <- function(df, hats)
{
  if (hats <  sum(df$MIN)) stop("Not enough hats to go round!")
  if (hats >= sum(df$MAX)) {df$HATS <- df$MAX; return(df)}
  df$HATS  <- df$MIN
  hats     <- hats - sum(df$MIN)
  while(hats)
  {
    weights  <- df$HATS/df$POPULATION
    allowed  <- which(df$HATS < df$MAX)
    smallest <- which.min(weights[allowed])[1]
    df$HATS[allowed][smallest] <- df$HATS[allowed][smallest] + 1
    hats <- hats - 1
  }

  return(df)
}

Now we try it with reasonable numbers:

dt1 %>% distribute_hats(50)
#>     ZONE POPULATION MIN MAX HATS
#>  1:  A34         40   1  10    3
#>  2: G345        110   0   9    8
#>  3:  H62         80   0   2    2
#>  4: D563         70   1  11    5
#>  5:  T63         90   0  12    6
#>  6: P983         90   1   8    6
#>  7:  S24        130   0   5    5
#>  8:  J54        140   1   3    3
#>  9: W953         80   1   2    2
#> 10:  L97         30   0   0    0
#> 11:  V56         80   1   8    6
#> 12:  R99         50   1   8    4

dt1 %>% distribute_hats(10)
#>     ZONE POPULATION MIN MAX HATS
#>  1:  A34         40   1  10    1
#>  2: G345        110   0   9    1
#>  3:  H62         80   0   2    1
#>  4: D563         70   1  11    1
#>  5:  T63         90   0  12    1
#>  6: P983         90   1   8    1
#>  7:  S24        130   0   5    0
#>  8:  J54        140   1   3    1
#>  9: W953         80   1   2    1
#> 10:  L97         30   0   0    0
#> 11:  V56         80   1   8    1
#> 12:  R99         50   1   8    1

And edge cases:

dt1 %>% distribute_hats(1000)
#>     ZONE POPULATION MIN MAX HATS
#>  1:  A34         40   1  10   10
#>  2: G345        110   0   9    9
#>  3:  H62         80   0   2    2
#>  4: D563         70   1  11   11
#>  5:  T63         90   0  12   12
#>  6: P983         90   1   8    8
#>  7:  S24        130   0   5    5
#>  8:  J54        140   1   3    3
#>  9: W953         80   1   2    2
#> 10:  L97         30   0   0    0
#> 11:  V56         80   1   8    8
#> 12:  R99         50   1   8    8

^{Created on 2020-05-09 by the reprex package (v0.3.0)}

I'm not sure about this. It sounds like an optimization or linear programming task

Here's the function:

allocate <- function(dt, N){
  if(N>dt[,sum(MAX)])
    stop("Too many hats to go around")

  if(N<dt[,sum(MIN)])
    stop("Not enough hats to go around")

# Allocate hats initially based on proportion but use clamping
  dt[, HATS := pmax(MIN, pmin(MAX, round(N * +(MAX>0) * POPULATION / sum(POPULATION[MAX>0]))))]

  n <- N - dt[,sum(HATS)]      
  if(n==0)  # All hats accouted for
    return(dt)

  if(n>0){  # Allocate the extra hats, again proportional to pop with clamping
    dt[HATS<MAX, HATS := HATS + pmax(MIN, pmin(MAX, 
              round(n * +(MAX>0) * POPULATION / sum(POPULATION[MAX>0]))))]
  } else {  # Or subtract the superfluous hats, according to pop
    dt[HATS>MIN, HATS := HATS - pmax(MIN, pmin(MAX, 
              round(abs(n) * +(MAX>0) * POPULATION / sum(POPULATION[MAX>0]))))]
  }

  n <- N - dt[,sum(HATS)]  # Check again
  if(n==0)  # All hats accouted for
    return(dt)

  if(n>0){  # This time, just add 1 hat to those that require them
    dt[HATS<MAX, i:=.I][i<=n, HATS := HATS + 1]
  } else {  # Or reduce the number of hats by one
    dt[HATS>MIN, i:=.I][i<=abs(n), HATS := HATS - 1]
  }

  dt[, i:=NULL]  # Remove this guy
  return(dt)
}

Test it for 50:

dt2 <- allocate(dt1, 50)
dt2
    ZONE POPULATION MIN MAX HATS
 1:  A34         40   1  10    2
 2: G345        110   0   9    8
 3:  H62         80   0   2    2
 4: D563         70   1  11    5
 5:  T63         90   0  12    7
 6: P983         90   1   8    7
 7:  S24        130   0   5    5
 8:  J54        140   1   3    3
 9: W953         80   1   2    2
10:  L97         30   0   0    0
11:  V56         80   1   8    5
12:  R99         50   1   8    4

50 hats were allocated.

It may not be elegant or mathematically sound, but that was my attempt for what it's worth. Hope it can be of some use.

Formulating this as an integer programming where objective function minimizing sum of squares between assigned and target allocation subject to minimum and maximum allocation constraints:

dt1[, TARGET := POPULATION / sum(POPULATION) * TOTAL]

system.time({
    library(CVXR)
    x <- Variable(nrow(dt1), integer=TRUE)
    mini <- dt1$MIN
    maxi <- dt1$MAX
    target <- dt1$TARGET
    obj <- Minimize(sum_squares(x - target))
    constr <- list(mini <= x, x <= maxi, sum(x) == TOTAL)
    prob <- Problem(obj, constr)
    result <- solve(prob)
})
#   user  system elapsed 
#   1.60    0.17    1.76 


dt1[, ALLOCATION := as.integer(round(result$getValue(x)))]

output:

    ZONE POPULATION MIN MAX   TARGET ALLOCATION
 1:  A34         40   1  10 2.020202          4
 2: G345        110   0   9 5.555556          7
 3:  H62         80   0   2 4.040404          2
 4: D563         70   1  11 3.535354          5
 5:  T63         90   0  12 4.545455          6
 6: P983         90   1   8 4.545455          6
 7:  S24        130   0   5 6.565657          5
 8:  J54        140   1   3 7.070707          3
 9: W953         80   1   2 4.040404          2
10:  L97         30   0   0 1.515152          0
11:  V56         80   1   8 4.040404          6
12:  R99         50   1   8 2.525253          4

data:

library(data.table)
dt1 <- data.table(ZONE = c("A34","G345","H62","D563","T63","P983","S24","J54","W953","L97","V56","R99"), 
  POPULATION = c(40,110,80,70,90,90,130,140,80,30,80,50), 
  MIN = c(1,0,0,1,0,1,0,1,1,0,1,1), 
  MAX = c(10,9,2,11,12,8,5,3,2,0,8,8))
TOTAL <- 50L