You can use limma::avereps (https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/avereps). To use this, you define a set of groupings over which to summarise (here each group would be contain a single Ensembl ID). However, avereps returns a matrix, rather than an ExpressionSet, so it's use would lead to you losing the phenotype / treatment data that you need in statistical modelling:

Here's my solution:

collapse_by_feature_ids <- function(eset, id = "Ensembl") {
  .filter <- function() {
    id_set <- Biobase::fData(eset)[[id]]
    keep_rows <- which(!is.na(id_set) & id_set != "")
    eset[keep_rows, ]
  }
  .collapse <- function(eset) {
    fdata <- Biobase::fdata(eset)
    gene_groups <- fdata[[id]]
    averaged_matrix <- limma::avereps(eset, ID = gene_groups)
    Biobase::ExpressionSet(                                                                        
      averaged_matrix,                                                                             
      phenoData = Biobase::phenoData(eset)                                                                  
    )
  }
  stopifnot(is(eset, "ExpressionSet"))
  stopifnot(id %in% colnames(Biobase::fData(eset)))
  filtered_eset <- .filter()
  .collapse(filtered_eset)                                                                                   
}

This ensures that the phenotype / treatment data is retained and that an expressionSet is returned. It takes the arithmetic mean of the expression-log-intensities for each Ensembl ID. You should add a feature data data-frame in afterwards