The multivariate Grey model (GM(1,N)) is a widely utilized forecasting tool within Grey models, capable of incorporating multiple variables. However, this model may present limitations in accuracy due to the structural differences between the Grey differential equations employed for parameter estimation and the whitening equation used for prediction. To overcome these issues, we propose an improvement to the GM(1,N) model by adjusting the residuals with a multivariate Markov chain. Thus, we develop the GMCM(1,N) model, which integrates the Mixture Transition Distribution (MTD) approach to model the interdependencies among residuals in multiple series. The use of the MTD model allows us to reduce the errors in parameter estimation in the Markov multivariate model. The proposed method improves the accuracy of representing dynamic interactions among variables and significantly increases the overall accuracy of the prediction. Finally, we apply the GMCM(1,N) model to the energy sector, a domain where multivariate Grey-Markov approaches are still scarcely explored. We test the performance of the model using historical data for both renewable and non-renewable energy production for Italy. Furthermore, accurate energy production forecasts could support strategic energy planning and policy development, highlighting the practical relevance and potential applications of the improved model within the Italian energy sector.