Zonal statistics for multiple rasters with overlapping polygons

I have a shapefile with 300 polygons that correspond to 300 rasters (the rasters have the exact same extent than their corresponding polygons).
What I’m trying to do is to run Zonal Statistics on this dataset but the problem is that many polygons overlap many rasters (the rasters can overlap each other as well).
Here is an example:

First I used this answer with the code below:

root = QgsProject.instance().layerTreeRoot()
for layer in root.children():
  if layer.name().startswith('RASTER_CODE'):
    prefix = layer.name()[-6:-4]
    params = {'INPUT_RASTER': layer.name(), 'RASTER_BAND': 1, 'INPUT_VECTOR': 'MY_SHAPEFILE', 'COLUMN_PREFIX': prefix+'_', 'STATS': 2}
    processing.run("qgis:zonalstatistics", params)

But since polygon 1 and polygon 2 overlap: the Zonal Statistic tool calculates 2 means for those polygons…

… and even if that part of the code:

if layer.name().startswith('RASTER_CODE'):
    prefix = layer.name()[-6:-4]

enables me to retrieve which "mean" is the correct one, it is far too time-consuming for 300 polygons…

Here what I would like to get:

So I could eventually sum up mean_1 mean_2 mean_3 in a new field mean gathering the mean of each of my polygons

EDIT: I found a workaround with Rstudio (see script below) but I would like to find a solution within pyQGIS

library(tidyverse)

data_with_row_names <- data.frame(data[,-1], row.names=data[,1])
means <- sapply(row.names(data_with_row_names), function(x) data_with_row_names[[x, x]])
means_final <- data.frame(t(sapply(means,c)))