Forecasting future heat load in smart district heating networks is a key problem for utility companies that need such predictions for optimizing their operational activities. From the statistical learning viewpoint, this problem is also very interesting because it requires to integrate multiple time series about weather and social factors into a dynamical model, and to generate models able to explain the relationships between weather/social factors and heat load. Typical questions in this context are: “Which variables are more informative for the prediction?” and “Do variables have different influence in different contexts (e.g., time instant or situations)?” We propose a methodology for generating simple and interpretable models for heat load forecasting, then we apply this methodology to a real dataset, and, finally, provide new insight about this application domain. The methodology merges multi-equation multivariate linear regression and forward variable selection. We generate a (sparse) equation for each pair day-of-the-week/hour-of-the-day (for instance, one equation concerns predictions of Monday at 0.00, another predictions of Monday at 1.00, and so on). These equations are simple to explain because they locally approximate the prediction problem in specific times of day/week. Variable selection is a key contribution of this work. It provides a reduction of the prediction error of 2.4% and a decrease of the number of parameters of 49.8% compared to state-of-the-art models. Interestingly, different variables are selected in different equations (i.e., times of the day/week), showing that weather and social factors, and autoregressive variables with different delays, differently influence heat predictions in different times of the day/week.