R ompr MILPModel array multiplication?

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() %>%
#total binaries ==10

#sum of binaries * mat1 <= 7
ompr::solve_model(with_ROI(solver = "glpk",verbosity=-2,gap_limit=0,time_limit=180, node_limit=-1,first_feasible=FALSE))$$`$$