Оптимизируйте среднее и стандартное отклонение групп в кадре данных, находя лучшее расположение

Question

Контекст

Я должен выделить 3 обработки с 12 повторениями для 36 участков в поле.Каждому участку присваивается переменная окружения, которая является фиксированной.

Задача

Как найти наилучшее возможное распределение моих обработок, чтобы для каждой обработки среднее и стандартное отклонение переменной среды было как можно более равным?

Данные

Некоторые фиктивные данные:

 df1 <- data.frame(plot= seq(1, 36, by=1),
              treatment =rep(1:3),
              enviroment=rnorm(36, 9300, 3000))
head(df1)
    plot var environment
1    1   1    8594.994
2    2   2    9737.577
3    3   3    8111.464
4    4   1    8367.244
5    5   2    5965.046
6    6   3    7664.372

Мой наивный подход

Это в основном просто грубое принуждение, но моя проблема, я полагаю, является комбинаторной и тяжелой для этого.подход.Чтобы подвести итог, я хочу минимизировать средние значения и стандартные отклонения для лечения.Посмотрите на Boxplots для уточнения.Я надеюсь, что это хотя бы немного интересная проблема, и буду признателен за любые идеи, как ее решить!Спасибо.

library(dplyr)    
# Boxplot by var to show distribution before optimization
    mean  <- df1 %>% group_by(var) %>% summarise(mean=mean(environment))
    sd    <- df1 %>% group_by(var) %>% summarise(sd=sd(environment)) %>% summarise(sum_sd=sum(sd))
    boxplot(environment~var,data=df1, main=paste("Before Optimization Sum-SD:",round(sd),"Diff_means:",round(abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2])))), 
            xlab="treatment", ylab="environment")

    # Run the randomization process n-times
    for(i in 1:10000000) {
      # Randomize treatment position over experimnetal units
      df1[,c("var")] <- df1[sample(nrow(df1)),c("var")]
      # Calculate Mean for each treatment
      mean <- df1 %>% group_by(var) %>% summarise(mean=mean(environment))
      # Calculate Sum of SD for all machines
      sd    <- df1 %>% group_by(var) %>% summarise(sd=sd(environment)) %>% summarise(sum_sd=sum(sd))
      # First run,save initial values
      if (i == 1) {
        # Save the difference of the means
        var_1 <- abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2]))
        # Save summed standard deviations
        sd_1  <- df1 %>% group_by(var) %>% summarise(sd=sd(environment)) %>% summarise(sum_sd=sum(sd))
      }
      # if randomized values are better then initial values save them
      if(var_1 > abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2]))&
         (sd_1 >sd)) {
        # save better values
        var_1 <-  abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2]))
        sd_1  <- sd
        # save improved experimental design
        best_setup <- df1
        # give visual output of better design
        boxplot(environment~var,data=best_setup, main=paste("nth Run:",i,"Sum-SD:",round(sd),"Diff_means:",round(abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2])))), 
                xlab="var", ylab="environment")
      }
    }
    # give visual output of best found experimental design 
    boxplot(environment~var,data=best_setup, main=paste("After Optimization Sum-SD:",round(sd),"Diff_means:",round(abs(abs(mean[1,2] - mean[2,2])+abs(mean[1,2] - mean[3,2])))), 
            xlab="var", ylab="environment")