R ompr MILPModel array multiplication?

Question

In R, I regularly ompr::MILPModel for optimization. I adapt the below snippet to enable multiplication of a decision variable with two dimensions (e.g., x[i,j] ) by a numeric matrix of the same dimensions, in constraints and the objective function. I chanced upon this code elsewhere, so I am not going to claim I much know what is going on matrix_multiplication_fcn, just that it works.

I would like to be able to use MILPModel with decision variables of 3+ dimensions, e.g., x[i,j,k] or x[i,j,k,m] , and be able to multiply these decision variables against numeric arrays of the same dimension. I am having a really hard time figuring out how to make an array_multiplication_fcn of 3+ dimensions that has the same effect.

I've made a couple attempts but when I look at the model's objective function, it's just the first two or three values from the numeric array, repeated over and over.

mat1 <- matrix(ncol=10,nrow=4,runif(400))

#define this function, it will be necessary for matrix multiplication inside a MILPModel
matrix_multiplication_fcn <- function(static_matrix, row_variable, column_variable){
  vapply(seq_along(row_variable), function(k) static_matrix[row_variable[k], column_variable[k]], numeric(1L))  }

milp_model <- ompr::MILPModel() %>%
  add_variable(assign_units[rowindex,colindex], rowindex=1:4,colindex=1:10,type='binary') %>%
  #total binaries ==10
  add_constraint(sum_expr( assign_units[rowindex,colindex],rowindex=1:4,colindex=1:10 )==10 ) %>%
  
  #sum of binaries * mat1 <= 7
  add_constraint( sum_expr( ompr::colwise(
    matrix_multiplication_fcn(static_matrix=mat1,row_variable=rowindex,column_variable=colindex)) *
                              assign_units[rowindex,colindex], 
                            rowindex = 1:4, colindex = 1:10) <= 7) %>%

  #objective: maximize value
  set_objective(sum_expr(
    ompr::colwise(matrix_multiplication_fcn(static_matrix=mat1,row_variable=rowindex,column_variable=colindex)) * 
      assign_units[rowindex,colindex],
    rowindex=1:4,colindex=  1:10),sense='max')
    
milp_model_out <-  milp_model %>%   
  ompr::solve_model(with_ROI(solver = "glpk",verbosity=-2,gap_limit=0,time_limit=180, node_limit=-1,first_feasible=FALSE))```