A large fraction of district heating networks in Europe are still operating without thermal energy storage. While offering operational flexibility, the presence of energy storage complicates the problems of optimal sizing and control of the heat production plant. This necessitates a novel approach to system sizing based on non-linear programming in order to capture the system’s thermal response and inertia. Therefore, the present study focuses on such a typical plant in France and describes a holistic approach to optimize the capacities/parameters of system components, positions of the energy sources, as well as hourly flow rates. Such a simultaneous tuning of the system parameters and operating variables is essential to maximize the cost savings over the system's lifetime. The proposed framework here maintains the original non-linear form of the objective function and the constraints for more accuracy while allowing precise modeling of the thermally stratified storage tank by decoupling its calculations at a lower simulation level. Hence, the computation time in MATLAB® R2020A is mostly below 20 min on a typical personal computer. The lowest levelized cost of heat (57.259 EUR/MWh) is achieved when placing the gas boiler on a parallel line between the biomass boiler and the storage tank while sizing the two boilers using fractions of 0.33 and 0.73 of the maximum hourly demand, respectively, as well as a storage volume of 318.9 m3. Such relatively small storage reduces the lifecycle costs by a maximum of 1.34 % while reducing the specific carbon emissions by 17.33 %. Besides, it dampens the fluctuations in supply flow rates, minimizes unnecessary flow recirculation, and ensures a stable heat supply. Overall, thermal storage has a reasonable positive impact from an economic perspective, but a substantial impact on operational flexibility and mitigated carbon emissions.