Linear Regression on data with bimodal outcome

I have a data set with 3,000 features and continuous dependent variables of time with 18,000 instances. The histogram of the dependent variables show that the they have a bimodal distribution. I am building linear regression models that forecast the time, but none of the models are able to make predictions; the $$R^2$$ values of all of the models are $$0$$. I plotted the residuals of the models and verified that they are normally distributed

I have used LassoCV/Lasso, ElasticNetCV/ElasticNet, XGBoost, and Kernel Ridge Regression in Scikit-Learn to make the models. I used the weights of the Lasso coefficients to filter my data frame so that I use the features that have the most magnitude in my predictions, but that did not work. I used the feature_importance_ method for XGBoost to do something similar, but that also did not help.

Kernel Ridge Regression ran with the parameters

KernelRidge(alpha=1, coef0=1, degree=3, gamma=0.001, kernel='polynomial', kernel_params=None)

but overfit the training data completely with $$R^2 = 1$$.

I have been looking into non-parametric methods like Kernel Density Estimate and Gaussian Mixture Models, but am uncertain in how to use them for my regression task. I have also been looking into Modal Regression and fitting of multi-functions, but am not familiar in how to implement them with Python.

I have also considered finding out what features are contributing to the modes and splitting the data set as two distributions and running the models on the two separate data sets, but am not sure on how to paint the picture for the entire data set.

In general, how should I approach building a regression model with a multimodal dependent variable?

One option would be use sklearn.compose.TransformedTargetRegressor to make the dependent variable more normal distributed.