The aim of this document is to present the topic of operational optimization in District Heating (OH) systems, with special focus on different kinds of thermal energy storage. An optimization solution based on solving multiple Mixed Integer Linear Programming (MILP) problems has been proposed and implemented in the R programming environment. The operational optimization in a DH system, especially if this system is supplied from a combined heat and power (CHP) plant, is a difficult and complicated task. Finding a global financial optimum requires considering long periods of time and including thermal energy storage possibilities into consideration. There are three important solutions for thermal energy storage: hot water tanks, utilization of thermal inertia of the network itself and utilization of thermal inertia of buildings. Each of these solutions has its advantages and disadvantages, and they can be combined to reach the maximum flexibility at lowest cost. However, modeling of operation with all of the thermal energy storage possibilities in place is a complicated task, since they influence the transient behavior of the network in different ways, and affect each other. On the other hand, optimal planning of heat production can be done only if simple and robust simulation models are available. Proposed solution allows simulation of three kinds of thermal energy storage, with their specific transient behaviors and interactions, at the same time keeping the model simple and ready to be used with a MILP solver. An iterative approach has been applied to non-linear phenomena, which allows solving a non-linear problem by multiple MILP optimizations. It has been successfully implemented in the "R" programming environment and tested on a simple example. The results can prove useful for DH system operators in the near future. (C) 2018 Elsevier Ltd. All rights reserved.