Operation hours optimization for circular schedule

Question

Here is my problem. A store has X = 15 electical devices with the ability to work non-stop, fully charged, up to 8 hours. Their battery charge lasts 2 hours and the operating hours of the store differ day over day, fortunately, this does not affect their charging. The number of chargers, Y = 4, most likely will be less than devices number. When the store remains closed charging takes place normally. Every device that is fully charged, remains charged if does not operate.Store remains open/closed according to operating hours (loo at the end)

Is their any way to formulate this problem into LP/MILP so as to achieve device operation hours maximization, for the indicated circular schedule?

Update - 7/6/2020

Following to @Richard' comment I trie to reduce the post and also share data to make possible any attempt for replication. I used R to formulate and get a first iteration.And ompr package to solve this.

Here is the code:

library(data.table)
library(stringr)
Test = data.table(permutations(n = 2,r = 12,v = c("b","s"),repeats.allowed = T))
Test[, names(Test) := lapply(.SD, function(x) gsub("b", "0,0", x))]
Test[, names(Test) := lapply(.SD, function(x) gsub("s", "1,1", x))]
Test[, key := do.call(paste, c(.SD, sep = ",")), .SDcols = names(Test)]
Conf = lapply(unique(Test$key),FUN = function(x) as.integer(str_split(x,",")[[1]]))

ConfDT = as.data.table(Conf)
ConfDT[, Hour := seq_len(.N)]
ConfDTMelt = melt(ConfDT,id.vars = "Hour",variable.name = "k",variable.factor = F,value.factor = F)
ConfDT[, Hour := NULL]
ConfDTMelt[, Change := rleid(value), .(k)]
ConfDTMelt[, CumValue := cumsum(value), .(k,Change)]
ConfDTMelt[, ValueLead :=shift(value,n = 1,"lead",fill = NA_real_), by = .(k)]
RmvCases = unique(c(ConfDTMelt[(Change > 1 & value == 1 & ValueLead == 0 & CumValue<8 & CumValue>0) | CumValue>8,k],ConfDTMelt[, .(Sum = sum(value)), by = .(k)][Sum == 0, k]))
ConfDT = ConfDT[, setdiff(names(ConfDT),RmvCases), with = F]
colnames(ConfDT) = paste0("V",seq_len(ncol(ConfDT)))
ConfDT[, Hour := seq_len(.N)]
ConfDTMelt = melt(ConfDT,id.vars = "Hour",variable.name = "k",value.name = "Work")
ConfDTMelt[, k:=as.integer(str_remove(k,"V"))]
ConfDT[, Hour := NULL]
Conf = as.list(ConfDT)

ConfDTComb = data.table(rbind(as.matrix(ConfDT[, c(CJ(1:ncol(ConfDT),1:ncol(ConfDT))$V1), with = F]),
                              as.matrix(ConfDT[, c(CJ(1:ncol(ConfDT),1:ncol(ConfDT))$V2), with = F])))
colnames(ConfDTComb) = paste0("Col",c(1:ncol(ConfDTComb)))
ConfDTComb[, Row := seq_len(.N)]
ConfDTCombMelt = melt(ConfDTComb,id.vars = "Row",variable.name = "Combination",value.name = "value")
ConfDTCombMelt[, Change := rleid(value), .(Combination)]
ConfDTCombMelt[, CumValue := cumsum(value), .(Combination,Change)]
ConfDTCombMelt[, ValueLead :=shift(value,n = 1,"lead",fill = NA_real_), by = .(Combination)]
ForbiddenCases = CJ(1:ncol(ConfDT),1:ncol(ConfDT))[as.numeric(str_remove(unique(ConfDTCombMelt[CumValue>8 | (CumValue<8 & ValueLead == 0 & CumValue >0 & Change>1),Combination]),"Col")),]

In the beginning I tried to get all possible permutations of c(0,1), which actually denotes when a devices is working or not. In order to reduce the total number of permutations I get the permutations for pairs of 0's and 1's. Finally I kept only the meaningful ones, so cases such as totally disabled/enabled devices or consecutive 0,1 where removed and I ended up with 105 configurations.

Model:

Fun1 = function(j,k){mapply(function(x,y) sum(unlist(Conf[y])*unlist(Rest[x])),j, k)}
ForbiddenCasesFun = function(k,kk){mapply(function(x,y) nrow(ForbiddenCases[V1 == x & V2 == y])==1 ,k, kk)}
ChargingHoursFun = function(k,h){mapply(function(x,y) unlist(Conf[x])[y] == 1,k, h)}

model = MILPModel() %>%
  # device- Day 
  add_variable(x[i,j], i = 1:15 , j = 1:7, type = "binary") %>%
  # device- Day - Configuration
  add_variable(y[i,j,k], i = 1:15, j = 1:7, k = 1:length(Conf), type = "binary") %>%
  # device- Day - Configuration - Hour charging
  add_variable(z[i,j,k,h], i = 1:15, j = 1:7, k = 1:length(Conf), h = 1:24,type = "binary") %>%
  # Single Configuration per device/Day
  add_constraint(sum_expr(y[i,j,k], k = 1:length(Conf)) <= x[i,j], i = 1:15, j = 1:7) %>% 
  # Day over Day forbidden adjacent configurations
  add_constraint(y[i,j-1,k] + y[i,j,kk] <= 1, i = 1:15, j = 2:7, k = 1:ncol(ConfDT),kk =  1:ncol(ConfDT),ForbiddenCasesFun(k,kk)) %>% 
  # Week over week forbidden adjacent configurations
  add_constraint(y[i,7,k] + y[i,1,kk] <= 1, i = 1:15, k = 1:ncol(ConfDT),kk =  1:ncol(ConfDT),ForbiddenCasesFun(k,kk)) %>% 
  # Charging can not be placed when devices operates according to the configuration  
  add_constraint(z[i,j,k,h] == 0, i = 1:15,j = 1:7, k = 1:ncol(ConfDT), h = 1:24, ChargingHoursFun(k,h)) %>%
  # single charging hour per configuration
  add_constraint(sum_expr(z[i,j,k,h],h = 1:24) == y[i,j,k], i = 1:15, j = 1:7, k = 1:length(Conf)) %>% 
  # Maximum number of devicescharging per hour
  add_constraint(sum_expr(z[i,j,k,h],i = 1:15, j = 1:7, k = 1:length(Conf)) <= 4, h = 1:24)
model = model %>%
  set_objective(sum_expr(y[i,j,k]*colwise(Fun1(j,k)),  i = 1:15, j = 1:7, k = 1:length(Conf))) %>%
  solve_model(with_ROI(solver = "symphony", control = list(verbosity = -1,time_limit = 500)))

Objective:

$$\sum_{\max} x_{i,j,k}\cdot \sum_{j,k} Conf_{k}\cdot Restr_{j} \ \ \ \forall i,j,k$$

Data Info

Operating Hours:

    Hour Mon Tue Wed Thu Fri Sat Sun
 1:    1   1   1   1   1   1   1   1
 2:    2   1   1   1   1   1   1   1
 3:    3   1   1   1   1   1   1   1
 4:    4   1   1   1   1   1   1   1
 5:    5   1   1   1   1   1   1   1
 6:    6   0   0   0   0   0   0   1
 7:    7   0   0   0   0   0   0   1
 8:    8   0   0   0   0   0   0   1
 9:    9   0   0   0   0   0   0   1
10:   10   0   0   0   0   0   0   1
11:   11   0   0   0   0   0   0   1
12:   12   1   1   1   1   0   0   1
13:   13   1   1   1   1   0   0   1
14:   14   1   1   1   1   0   0   1
15:   15   1   1   1   1   0   0   1
16:   16   1   1   1   1   0   0   1
17:   17   1   1   1   1   0   0   1
18:   18   1   1   1   1   0   0   1
19:   19   1   1   1   1   0   0   1
20:   20   1   1   1   1   0   0   1
21:   21   1   1   1   1   0   0   1
22:   22   1   1   1   1   0   0   1
23:   23   1   1   1   1   1   1   1
24:   24   1   1   1   1   1   1   1

What I see is that the problem is really heavy, even model building is highly time consuming. What I have not figure out yet is how is that even possible to define a configuration-based charging. If configuration starts with 1's then the charging should take place in the previous day's configuration. And if we set different charging properties due to the load consumption how can this be formulated. Charging may last longer in the noon where light and A/C are on and capacity remains the same, then charging may last 2-3 hours. Any suggestions?