I'm a 22 year old Statistics student with a big problem to solve. English isn't my first language, so I apologize in advance for any mistakes in my grammar.
I'm trying to make a short-term gas consumption forecast for a small city. In my country, gas consumption has its peak in the winter months (may, june, july, august and september), and in the rest of the year the consumption is much lower.
I'm working with daily data from 2015 up until 2022, but I can, if necessary, disregard the first few years. Since I have daily data (one observation of the gas consumption per day), I have weekly and yearly seasonality. I also need to take into account some explanatory variables, the most important one being the daily mean temperature of the city.
Here I include a picture of the time series across the years.

To sum it up, I need a forecasting model that will handle daily data with two types of seasonalities (weekly and yearly) and that also allows me to incorporate covariates (the temperature and also some dummies for holidays and weekends). In this context, with what type of model should I work with?
After applying a Box-Cox transformation to stabilize the variance of my variable (gas consumption), I tried using dynamic regression, following Rob Hyndman's book. I incorporated the following variables as covariates in the model: a dummy for saturdays, another for sundays and holidays, the mean daily temperature (in its linear and quadratic form) and the heating degree days (in its linear and quadratic form). The residuals plots look like this:

I also tried a dynamic harmonic regression with the following covariates: fourier terms, the mean daily temperature (in its linear and quadratic form) and the heating degree days (in its linear and quadratic form). The code I used was:
fit7 <- time_series %>%
model(
ARIMA( gasconsumption_transformed ~ PDQ(0, 0, 0) + pdq(d = 0) +
mean_temp + I(mean_temp^2) + hdd + I(hdd^2)+
fourier(period = "week", K = 3)+
fourier(period = "year", K = 5))
)
report(fit7)
fit7 %>% gg_tsresiduals()
And the residual plots were oddly similar to the previous fit with dynamic regression.

The residuals of the two previous models don't look particularly good. So I'm hoping to find another possible solution to my problem.
I've also read Rob J. Hyndman's paper titled "Short-term load forecasting based on a semi-parametric additive model", but I haven't found any package that may help me do this in the software R.
Another type of model that I've come across is state space models however in this case I haven't found any R package that implements this method either.
Any type of help with models that are easily implemented in R would be greatly appreciated! Also, if there's something hard to understand or information missing, I'll gladly re-check the cross-validated question.
One final comment: I am not allowed to share the data used.