En este trabajo, derivamos la matriz operacional utilizando polinomios de Bernoulli poliédricos. Estos polinomios generalizan los polinomios de Bernoulli utilizando una función generadora que involucra una función de polilogaritmo. Primero mostramos algunas propiedades nuevas para estos polinomios de Bernoulli poliédricos; luego derivamos una nueva matriz operacional basada en polinomios de Bernoulli poliédricos para la derivada de Atangana-Baleanu. Se deriva una matriz operacional de retardo basada en polinomios de Bernoulli poliédricos. Se muestra el límite de error de este nuevo método. Aplicamos esta matriz operacional de Bernoulli poliédrica para resolver ecuaciones diferenciales fraccionarias de retardo con coeficientes variables. Los ejemplos numéricos muestran que este método es fácil de usar y aún así capaz de dar resultados precisos.