The aim of this study is to present an optimal demand response (DR) program using the IEC 60076-7 standard thermal model in order to simultaneously determine the effect of operating conditions and loading on the thermal status and failure rate of power transformers, as well as to determine the optimal point for changing the operating mode (single or parallel) based on economic criteria including the Expected Energy Not Supplied (EENS) and power losses. None of the previous references have considered basic reliability indices (such as the failure rate) as the main criterion for decision-making regarding transformer replacement. Furthermore, the effect of loading on the failure rate, as well as the effect of specific operating conditions in which the hot-spot temperature exceeds the permissible limit, have not been taken into account in previous models.Presenting a new mathematical relationship to model the effect of loading on the transformer failure rate, completing the mathematical model of the failure rate by incorporating the effect of loading on the hot-spot temperature, and formulating a linear optimization problem (utilizing piecewise linearization of nonlinear thermal functions) that simultaneously optimizes the volume and management of demand response over a specified time horizon while considering thermal constraints (hot-spot and oil temperatures), maximum current (1.5 per unit), and insulation aging (AEQ ≤ 1) are among the most important innovations of this paper.The achievable reserve margin through applying DR is up to 240 kVA, which is equivalent to 3.4 times the traditional reserve (70 kVA), without violating the thermal limit of 120 °C. In the "energy shift" mode, only 64% of the estimated DR capacity is utilized. The DR volume can reach up to 30% of the peak load, but this is only required under N-1 conditions and for a few hours per year. Under normal (N) conditions, no demand response is required, and no thermal stress is observed for the transformer, even with a high reserve margin. However, the IEC 60076-7 standard thermal model, due to its high temporal resolution (1 minute), causes a dramatic increase in variables and constraints over long-term horizons, rendering the optimization problem unsolvable. Piecewise linearization of the nonlinear thermal functions significantly reduces computation time compared to nonlinear models, but it remains challenging for long time periods (several days to a year). The proposed method effectively balances loading and thermal constraints, freeing up reserve capacity that was previously inaccessible; however, further research is required for long-term horizons.