Should I use Linear Regression in this case? I'm doing a personal project on regression and I hope I can get some advice from you on several problems. The dataset I'm having is about cuisine, shape is 180x9, and the 2 continuous variables are "prep_time" and "cook_time". I want to predict cook_time based on other variables. I one hot encoded other categorical variables and get a total of ~30 features in the end (and only prep_time is continuous). This is the dataset before processing:

The baseline model with linear regression, SVR and Random forest doesn't perform really well (The absolute r2 scores are <0.1 for all 3 models). But when I added a feature total_time = prep_time+cook_time, and used MinMaxscale on prep_time and total_time, the RF model improves to r2 of ~0.3. However, LinReg results in ~0 mse and 1.0 r2 on the test set, surprisingly.
Thus, I tried to analyze this LR model more. I used sklearn.feature_selection.SelectFromModel function to get the important features from this Linear Regression model, and it shows that only the prep_time and total_time are important for the model. In fact, the model trained using only these 2 variables perform just as well. VIF of these 2 variables are ~7.0, which implies multicollinearity. However, if I remove either of these 2 variables, and train on only the other feature, the model doesn't perform that well.
My question is, is this still a good model? What could explain/potentially be a problem with this LinReg model when it has such a low mse and high r2? Should I use another algorithm instead?
Thank you!

I suspect the issue is:

But when I added a feature total_time = prep_time + cook_time

Your response variable is cook_time. So it is nonsense to include it as a predictor. Suppose your model is something like:
$$y = x_1 + x_2 + x_3 + x_4$$
where $y$ is cook time, and $x_1$ is prep time.
So let us introduce another variable to represent total time: $z = y + x_1$,
or $x_1 = z - y$.  Now we substitute this into the model and we get:
$$y = z - y + x_2 + x_3 + x_4$$
So now we have $y$ on both sides of the model equation. This is going to lead to a perfect fit, which means that $R^2$ will be 1.0) and there will be no residuals (so that $MSE$ will be zero).
Moreover, let's just say that you did proceed with such a model. Now, when it comes to making predictions on new (out of sample) data, for which you obviously don't know what the cook time is, how can you form the total_time feature to predict on the new data ?
