Convergence issues and model selection in glmmTMB Convergence problems in mixed effect models seem to be a common struggle. It is my understanding that they emerge when the likelihood surface is too flat for the optimisation algorithms to find a single set of parameters that maximises the likelihood. Three of the common approaches to address them is to rescale the explanatory variables, try different optimizers, and/or simplify the model (see here, here, and here). I have followed these suggestions but they still led me an impasse: namely the model that reflects the nature of the data best does not converge, while a simpler model converges but it does not capture the variability associated to the random factors as accurately.
What follows is a detailed reproducible example to illustrate my question.
Let's say I am looking at how unicorn herd size changes with food quantity. Unicorn herd sizes were surveyed at seven localities over the course of twelve months. Food quantity was assessed monthly, so for each locality at each month we have one value of food quantity (continuous predictor, let's say grass cover). For each locality on each month, unicorn herds were re-surveyed over four consecutive days to account for both unicorn movement and possible sampling error. On each sampling day I have a variable number of data points, each referring to the size of a different herd. Here is a diagram of the data structure (Fig.1):

Simulated data available here:
unicorns <- read.csv(text= RCurl::getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/unicorns.csv"))

Here is an overall look at the data (Fig.2):

Looking at the data separately by locality it appears that the relationship between food availability and herd size differ among localities (Fig.3):

The relationship between food availability and herd size differs also among survey months, but by eyeballing the graphs it looks like survey month is a smaller source of variability compared to locality (Fig.4):

Attempting to fit a model to the data with random locality (including slopes), random sampling month (including slopes), while also accounting for temporal autocorrelation of repeated measures fails:
unicorns$survey_day_f <- glmmTMB::numFactor(unicorns$survey_day_n)
levels(unicorns$survey_day_f)

unicorns_glmmTMB0 <- glmmTMB::glmmTMB(Herd_size_n ~ food.quantity
                        + (1 + food.quantity | Locality)
                        + (1 + food.quantity | Year_Month)
                        + ou(survey_day_f + 0 | Locality),
                        family="poisson",
                        data=unicorns)
# Error in eigen(h) : infinite or missing values in 'x'

Refitting the model after scaling food.quantity works but doesn't converge:
dd.CS <- transform(unicorns, food.quantity =scale(food.quantity))
unicorns_glmmTMB.cs <- glmmTMB::glmmTMB(Herd_size_n ~ food.quantity
                        + (1 + food.quantity | Locality)
                        + (1 + food.quantity | Year_Month)
                        + ou(survey_day_f + 0 | Locality),
                        family="poisson",
                        data= dd.CS)
# Warning messages:
# 1: In fitTMB(TMBStruc) :
  # Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
# 2: In fitTMB(TMBStruc) :
  # Model convergence problem; singular convergence (7). See vignette('troubleshooting')

Changing optimiser does not help:
unicorns_glmmTMB.csb <- update(unicorns_glmmTMB.cs, control=glmmTMBControl(optimizer=optim,
               optArgs=list(method="BFGS")))
# Warning message:
# In fitTMB(TMBStruc) :
#   Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

A simpler model, which allows for random slopes for localities but not for survey months, still does not converge:
unicorns_glmmTMB <- glmmTMB(Herd_size_n ~ food.quantity
                            + (1 + food.quantity | Locality)
                            + (1 | Year_Month)
                            + ou(survey_day_f + 0 | Locality),
                            family="poisson",
                            data= unicorns)

Rescaling, changing optimiser to BFGS, or doing both did not help.
The next step of model selection/simplification led me to a model allowing only for random intercept for both Locality and survey month (Year_Month):
unicorns_glmmTMB_noslope <- glmmTMB::glmmTMB(Herd_size_n ~ food.quantity
                        + (1 | Locality)
                        + (1 | Year_Month)
                        + ou(survey_day_f + 0 | Locality),
                        family="poisson",
                        data= unicorns)
unicorns_glmmTMB_noslope$sdr$pdHess ## there is convergence

This models converges and does not show anything dodgy. The problem is that, from figure 3, I can tell that a model including (1 + food.quantity | Locality) would better reflect the nature of my data compared to a model including (1 | Locality).
I am stuck between an accurate model that gives unreliable estimates and a less accurate one that gives reliable estimates. Is there a solution to this conundrum that I am not seeing?
 A: Another option, suggested by Robert LaBudde and outlined by Ben Bolker here, is to treat problematic random effects as fixed effects. In the case of Locality it makes sense because differences among localities might reflect ecological patterns or might result from ecological processes. Adding Locality as a fixed factor can be done in two ways: by allowing it to have an interactive effect with food.quantity (interactive model) or not (additive model).
# Interactive model:
unicorns_glmmTMB_fixedloc <- glmmTMB::glmmTMB(Herd_size_n ~ food.quantity*Locality
                        + (1 | Year_Month)
                        + ou(survey_day_f + 0 | Locality),
                        family="poisson",
                        data= unicorns)
unicorns_glmmTMB_fixedloc$sdr$pdHess ## there is convergence

# Additive model:
unicorns_glmmTMB_fixedloc_additive <- glmmTMB::glmmTMB(Herd_size_n ~ food.quantity+Locality
                        + (1 | Year_Month)
                        + ou(survey_day_f + 0 | Locality),
                        family="poisson",
                        data= unicorns)
unicorns_glmmTMB_fixedloc$sdr$pdHess ## there is convergence

The two models can then be compared (with caution):
AICc(unicorns_TMB_fixloc, unicorns_TMB_fixloc_additive)
#                              df     AICc
# unicorns_TMB_fixloc          17 5792.111
# unicorns_TMB_fixloc_additive 11 5794.318
anova(unicorns_TMB_fixloc, unicorns_TMB_fixloc_additive)
#                              Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
# unicorns_TMB_fixloc_additive 11 5793.8 5840.6 -2885.9   5771.8                           
# unicorns_TMB_fixloc          17 5790.9 5863.2 -2878.4   5756.9 14.905      6    0.02101 *

I think that the outcome of model selection between unicorns_glmmTMB_fixedloc and unicorns_glmmTMB_fixedloc_additive can be considered a proxy of what the model selection outcome would have been between unicorns_glmmTMB and unicorns_glmmTMB_noslope, had unicorns_glmmTMB converged. In this case we are still on the fence (cfr. the significance of the ${\chi}^2$ test versus the small $\Delta$AIC and $\Delta$AICc), but in a way that's reassuring: it suggests that by settling for unicorns_glmmTMB_noslope over unicorns_glmmTMB we are missing out on a small amount of information.
