**2. End-Users Side Flexibility for Grid Upgrades Deferral**

We suggest the reader to start reading from Section 2.1 presenting the case study of this paper. Moreover, this section explains the choice of assumptions allowing us to consider that the calculated DR has safety margins to mitigate the thermal constraints of transformers. Further, Section 2.2 is intended to help the reader understand at what reserve margins it would be sufficient to use DTR only and from what moment it is necessary to apply the DR and DTR. Finally, in Section 2.3, the reader can see an overview of the proposed methodology to find the necessary DR.

#### *2.1. Case Study*

The case study (Figure 2a) is an outdoor MV/LV substation with 2 × 500 kVA distribution transformers both equipped with an ONAN cooling system (Oil Natural Air Natural). Figure 2b shows the considered annual load profile at 1-h resolution representing an aggregated consumption of one hundred houses simulated by physical and behavior approaches [46].

To extract the load profile, authors used a MATLAB application "House load" [47] with all default parameters except the ambient temperature profile (*θ<sup>a</sup> <sup>t</sup>* ). The updated hourly *θ<sup>a</sup> <sup>t</sup>* in Grenoble, France was provided by MeteoBlue for 2019 [48]. To obtain the conservative reserve values, only the historical monthly maximum for the last 30 years was used in thermal simulations (Figure 2b) [23]. Historical maximums of *θ<sup>a</sup> <sup>t</sup>* obtained for each month are assumed as a safety margin when simulating the worst thermal state of the transformer for different reserve values. This is especially relevant for global warming with expected constantly rising ambient temperatures. Moreover, this margin can mitigate other errors in thermal characteristics/modeling.

An additional margin in the study comes from the computation of the reserve that assumes an N-1 condition with the loss of one transformer in the considered substation. Thus, if one transformer is out of service, then the remaining distribution transformer 500 kVA can definitely supply all load alone, since the registered peak power is only 430 kVA (86% of nominal rating, 500 kVA). Therefore, the reserve for load connection (in a traditional DSO approach) would be estimated as the difference between the nominal rating and the peak load i.e., 500 − 430 = 70 kVA (i.e., 14% of nominal rating). Instead of the nominal rating, the DSO can also apply other approaches as static thermal ratings [49,50], seasonal ratings [50], or emergency ratings [51]. Nevertheless, the approach for reserve determination does not change in its nature, as it still represents a simple difference between admissible constant rating and the peak load without consideration of a true thermal state of transformers.

Once thermal modeling is applied, the transformer temperatures (i.e., oil temperature *θo <sup>t</sup>* and hot spot temperature *θ<sup>h</sup> <sup>t</sup>* ) and the Aging EQuivalent (AEQ) over the year can be calculated explicitly using the loading profile *Ptr <sup>t</sup>* (in p.u.) and ambient temperature. Current, temperature, and aging, as any variable representing the physical state, have corresponding limits. The choice of those limits can be again used to set important safety margins for thermal modeling. For instance, Table 1 provides current and temperature limits for various loading types: normal cyclic loadings as well as two emergency modes long-term overloading and short-term overloading. Normal cyclic loading refers to a situation when a transformer is subject to high ambient temperature or higher-than-rated load, but the thermal aging remains the same as for nominal conditions. Long-term emergency loading is a situation when a transformer is subject to elevated temperatures for days or even months, and the short-term emergency loading is a heavy overloading of a transformer for less than 30 min. Due to the temporal nature of emergencies, the IEC standard [42] allows increasing their temperature and current limits.


**Table 1.** Limits for distribution transformers applicable for different types of loadings [42].

In this paper, we use the normal cyclic loading as the strictest limits in the N-1 condition. However, DSO could choose the long-term emergency limits as an alternative for N-1 conditions. This would allow releasing more transformer capacity at the cost of higher risks of overheating and accelerated loss of life. This alternative is not considered in the paper but could be easily integrated without changing the proposed methodology and algorithms.

As stated in [23], the problem of reserve determination is that the load profile of a new consumer cannot be definitely known in advance. Therefore, the transformer's load profile after the connection of new consumers is unknown as well. In addition to the renewable production, leading to the famous duck curve from California, it is believed that the electrification of the heating and transport sector will likely change typical shapes of existing load profiles. In order to add a further "safety margin", the existing shape of a load profile is not increased proportionally in this paper. Instead, the reserve is considered while adding a constant load to the existing load profile all along the representative year [23]. That constant added load profile ensures the worse thermal mode of operation in comparison to any other pattern with the same peak power. Nevertheless, this assumption considers the peak increase only from new load connections but existing consumers, especially industry customers, can also increase the peaks. DSO should still guarantee that consumers can withdraw all power from the distribution network, which is indicated in connection contracts (also known as firm capacity contracts). In fact, consumers usually do not use their full allocation of power, but there is no guarantee that one day, some industrial or commercial consumers will not decide to expand production capacities and then boost the actual power demand up to the subscribed power indicated in connection contracts (of course without informing DSO). In such a case, the above-mentioned assumption can lead to errors in substation peak estimations and even to emergencies. To mitigate such risks, DSO could add a full-subscribed power (i.e., not measured power) of large industrial and commercial consumers to the existing load profile.

#### *2.2. Problem Statement*

The problem statement is described in Figure 3 by conducting preliminary thermal studies. These preliminary studies investigate what would be *θh, θo,* and AEQ for different reserve margins without using DR (i.e., without taking any measures to decrease the temperatures of transformers). Figure 3 displays the state variables of the transformer as a function of the added constant load from 1% to 100% of a nominal rating of 500 kVA to the existing load profile shown in Figure 2b. The maximal *θh, θ<sup>o</sup>* during the year, and corresponding AEQ are estimated for the N-1 condition and computed using the IEC 60076-7 standard [42]. The IEC 60076-7 standard is an internationally recognized loading guide that provides mathematical equations and thermal characteristics of oil-immersed transformers needed to calculate *θh, θo,* and AEQ. Specifically, we use the difference equations described in the annex E of the IEC 60076-7 standard. These equations are presented later in Section 3.1 of this paper.

The preliminary thermal studies show that it is possible to connect a constant load of 240 kVA (=0.49 p.u \* 500 kVA) in addition to the existing consumption without violating the thermal constraints at 120 ◦C (see point 1). That 240 kVA reserve is 3.4 times higher than the reserve (70 kVA) calculated with a conventional approach—i.e., the transformer rated power (500 kVA) minus the peak consumption (430 kVA). The 240-kVA reserve is even more significant than the similar reserve obtained for DTR based on design temperature 98 ◦C (145 kVA, 0.29 p.u. at point 0). Although the use of design temperature in DTR is often claimed to avoid the accelerated aging, we observe that the accelerated ageing occurs only if reserve (a load growth) reaches significant values around 0.70 p.u. (point 3). Thus, the consideration of *θ<sup>h</sup> <sup>t</sup>* limit should be preferred over design temperature if the aging limit is explicitly taken into account [23].

**Figure 3.** Preliminary thermal results for yearly simulation: load growing from 1% to 100% of nominal rating (reserve in p.u.). The initial loading of the studied distribution transformer is 86% (see the *y*-axis of the middle figure).

Assuming that appropriate DR programs prevent the violation of temperature constraints, more load can be further connected to the transformer until the next limit is reached—current limit at 1.5 p.u. (point 2 at 350 kVA in Figure 3). Starting from this point, DSO should use DR to avoid breaking both *θ<sup>h</sup> <sup>t</sup>* and current limits. Thus, from point 2 onward, the reserve can be further increased from 320 to 350 kVA until the next critical point is reached (see point 3 in Figure 3). Starting from this point 3, the AEQ may be higher than the normal annual loss of transformer life. Hence, it is necessary to reduce the *θ<sup>h</sup> t* even well below its limit to keep the AEQ less than 1 (as it will be shown in Section 4). The last critical point 4 in Figure 3 (reserve of 365 kVA) is a situation when the transformer oil temperature violates its own limit. Finally, an appropriately designed and managed Demand Response should tackle three constraints at the same time—transformer thermal limits (*θ<sup>o</sup> <sup>t</sup>* , *θ<sup>h</sup> <sup>t</sup>* ), current limit, and insulation ageing—which is the objective of the optimization problem introduced in Section 3. Note that heavy computational burden may arise due to the needs for the thermal modeling in accordance with IEC 60076-7 that sets specific requirements for the time step of the load profiles and the temperatures. Indeed, IEC standard [42] requires that the time resolution must be at least two times smaller than the winding time constant, which is usually about 4–10 min. Otherwise, the thermal calculations may lose numerical stability. Practically, this requires that the time step of the thermal model for the transformer does not exceed 2–5 min or may be even less. In this paper, a 1-min resolution is adopted.

It is important to note that DR will be required only a few days per year (e.g., high load, high ambient temperature, or maintenance works). Thus, the DR design and management can be formulated deterministically only for the days when the transformer violates the thermal limits and not for the entire year. Hence, we avoid formulating a large intractable optimization problem. However, we argue that this idea remains valid if the longest interval does not exceed a few days. Figure 4 shows the longest interval with thermal violations considered in the DR optimization problem as a function of the reserve margin

(added constant load on the *x*-axis). From Figure 4, we see that the longest interval increases exponentially. For a reserve value of 1.04 p.u. onward, the overheating occurs every day of the simulated year. Hence, nonlinear equations from IEC 60076-7 would lead to the intractability of the optimization problem. Therefore, in Sections 3.3 and 3.4, the paper suggests a linearization of the nonlinear equations of the transformer thermal model.

**Figure 4.** Number of days where thermal limits are reached.

#### *2.3. Methodology*

Figure 5 shows the algorithm that allows us to compute the DR required to ensure a given level of reserve.

The workflow consists of three main stages. At the first stage, the annual transformer loading (without DR) is computed with the existing load profile increased by the given reserve margin. Furthermore, the thermal state of the transformer is estimated with the ambient temperature over the whole year. If there are no thermal violations, then no DR is needed, and the algorithm stops here. Hence, another reserve margin can be investigated.

If any overheating is detected, the second stage of the algorithm identifies the interval(s) where transformer temperatures (*θ<sup>o</sup> <sup>t</sup>* , *θ<sup>h</sup> <sup>t</sup>* ) or current limits are violated. The algorithm extracts the loading profile (*P<sup>l</sup> <sup>t</sup>* ) at every identified interval. Then, the load and ambient temperature profiles over a given interval are considered as inputs for the integrated DR management and design that computes the minimum DR needs to fulfill the operating constraints (see Section 3). Note that the thermal model of the transformer requires initial values for the top oil and hot-spot temperatures at the beginning of the extracted interval. To do that, the optimized transformer loading profile (i.e., with DR management) from previous calculation is used to update the annual load profile from the start of the year until the beginning of the next interval. This cycle repeats until all intervals are investigated and allows us to track the correct initial temperature every time an interval is simulated. The last stage of the algorithm is needed once all intervals are studied and the optimized annual load profile is entirely reconstructed. Finally, the algorithm defines the DR values in terms of power and energy ratings as the maximum DR power and energy computed over the different intervals.

**Figure 5.** The procedure for finding the needed Demand Response (DR) volume to interconnect the studied reserve.

#### **3. Integrated Design and Management of Demand Response**

To understand the thermal modeling of the chosen ONAN transformer, the reader can refer to Section 3.1. Section 3.2 suggests the problem formulation of integrated design and management for DR. As we mentioned earlier, the nonlinear thermal equations of IEC 60076-7 would lead to the intractability of the optimization problem at long intervals. That is why Sections 3.3 and 3.4 explain how nonlinear equations have been linearized in this paper.

#### *3.1. Transformer Thermal Model and Aging*

The thermal model of oil-immersed transformers as well as the values for thermal characteristics are derived from the IEC 60076-7 standard [42]. The IEC model allows a discrete representation of the differential equations that govern the thermal behavior of the transformer. The specific thermal characteristics of studied ONAN distribution transformers are given in Table 2. The meaning of each symbol used either in Table 2 or in the next equations is provided in Table 3. Note that the values of Table 2 are very conservative, since IEC 60076-7 was intended to represent the whole transformer fleet by using the same set of thermal characteristics [52]. Generally, the thermal characteristics of transformers are obtained by performing so-called non-truncated heat run tests [53]. Non-truncated heat run tests mean a situation when a constant load is applied to the transformer until reaching the steady-state (constant) temperatures of oil and winding [54]. The reader can see [52,53] for more details on thermal characteristics and the equations of IEC 60076-7 [42].


**Table 2.** Thermal characteristics of ONAN (Oil Natural Air Natural) distribution transformers at nominal conditions.

**Table 3.** List of the used symbols.


The IEC model estimates the top oil and hot-spot temperatures over time, *θ<sup>o</sup> <sup>t</sup>* and *θh <sup>t</sup>* , respectively. As already mentioned, those values depend on the time-series profiles for the ambient temperature (*θ<sup>a</sup> <sup>t</sup>* ) and the transformer loading (*Ptr <sup>t</sup>* in p.u.) as well as the transformer thermal characteristics (Table 2). Specifically, the model includes two nonlinear functions denoted f 1 *<sup>t</sup>* and f 2 *<sup>t</sup>* (1), which are used within equations for *θ<sup>o</sup> <sup>t</sup>* and *θ<sup>h</sup> t* . The equations of IEC standard are ultimately summarized in (2) (for *t* > 1 min) and in (3) (for initialization step i.e., *t* = 1 min).

$$\mathbf{f}\_{l}^{1}\left(P\_{l}^{tr}\right) = \boldsymbol{\Delta\theta}^{cr} \times \left(\frac{\mathbf{1} + \left(P\_{l}^{tr}\right)^{2} \times \boldsymbol{R}}{1 + \boldsymbol{R}}\right)^{\mathbf{x}} \text{ and } \mathbf{f}\_{l}^{2}\left(P\_{l}^{tr}\right) = \boldsymbol{\Delta\theta}^{lr} \times \left(P\_{l}^{tr}\right)^{\mathbf{y}}\tag{1}$$

$$\begin{cases} \begin{aligned} \theta\_{t}^{o} &= \theta\_{t-1}^{o} + \frac{1}{k^{11} \times \tau^{o}} \Big( \mathbf{f}^{1} \left( P\_{t}^{tr} \right) + \theta\_{t-1}^{o} - \theta\_{t}^{o} \Big) \\ \Delta \theta\_{t}^{b1} &= \Delta \theta\_{t-1}^{b1} + \frac{1}{k^{21} \times \tau^{w}} \times \left( k^{21} \times \mathbf{f}\_{t}^{2} \left( P\_{t}^{tr} \right) - \Delta \theta\_{t-1}^{b1} \right) \\ \Delta \theta\_{t}^{b2} &= \Delta \theta\_{t-1}^{b2} + \frac{k^{22}}{\overline{\tau}^{v}} \times \left( \left( k^{21} - 1 \right) \times \mathbf{f}\_{t}^{2} \left( P\_{t}^{tr} \right) - \Delta \theta\_{t-1}^{b2} \right) \\ \theta\_{t}^{l} = \theta\_{l}^{o} + \Delta \theta\_{t}^{vl1} + \Delta \theta\_{t}^{lv2} \end{aligned} \end{cases} \quad \forall t \in T$$

$$\begin{cases} \theta\_{t=1}^{o} = \dot{\mathbf{f}}^{2} \left( P\_{t=1}^{tr} \right) + \theta\_{t=1}^{a} \\ \Delta \theta\_{t=1}^{k1} = k^{21} \times \mathbf{f}\_{t=1}^{2} \left( P\_{t=1}^{tr} \right) \\ \Delta \theta\_{t=1}^{h2} = \left( k^{21} - 1 \right) \times \dot{\mathbf{f}}\_{t=1}^{2} \left( P\_{t=1}^{tr} \right) \end{cases} \tag{3}$$

The thermal model also returns the annual equivalent aging (denoted AEQ) of the transformer insulation for the given ambient temperature and power profiles. The insulation degradation is computed according to (4) with the hot-spot temperature value over the simulated horizon (i.e., *t*∈*t*) at 1-min resolution here. Note that the aging is normalized with the period duration (the cardinal function #*T*) and has to remain below 1, which corresponds to a normal degradation with the transformer operating at the design temperature along its estimated lifetime (see the threshold on AEQ violation in the algorithm in Section 2).

$$\text{AEQ} = \frac{1}{\#T} \times \sum\_{t \in T} 2^{\frac{\vartheta\_t^h - \vartheta \aleph}{6}} \le 1 \tag{4}$$

#### *3.2. Problem Formulation*

The characterization of the DR volume, which is necessary to avoid the transformer overheating, is expressed in the form of a systemic optimization problem. In this optimization problem, the management strategy of the DR (i.e., load power profile modification) is considered along with DR design (sizing). Then, both sizing and management are variables of a single problem. Moreover, dynamic constraints should be considered due to the time dependency of the temperature profiles. The overall problem simply consists of minimizing the DR needs in terms of rated power (*PDRr* in kW) and the rated capacity (*EDRr* in kWh). The DR should ultimately fulfill the transformer thermal, aging, and loading constraints given in (5). Additional constraints are introduced to represent the DR operation *PDR t* within its bounds (6). This management allows us to modify the transformer loading (*Ptr <sup>t</sup>* <sup>×</sup>*Ktr*) for a given load profile (*P<sup>l</sup> <sup>t</sup>* ) following the power balance constraint in (7).

$$obj: \min \left( P^{DRr} + \frac{E^{DRr}}{1h} \right) \\
\text{s.t.} \begin{cases} \theta\_t^o \le \overline{\theta^o} \quad \forall t \in T \\ \theta\_t^h \le \overline{\theta^h} \quad \forall t \in T \\ P\_t^{tr} \le 1.5 \text{ p.u} \ \forall t \in T \\ \underline{AEQ} \le 1 \end{cases} \tag{5}$$

$$t - P^{DRr} \le P\_t^{DR} \le P^{DRr} \quad \forall t \in T \tag{6}$$

$$P\_t^{tr} \times K^{trK} = P\_t^l - P\_t^{L\mathcal{H}} \quad \forall t \in T \tag{7}$$

So far, no economic criteria are taken into consideration and no cost is attached to the capacity of the DR flexibility (e.g., cost for storage capacity) and its power (e.g., cost for a battery inverter or backup generator). Note that in practice, this DR flexibility could be of any form and provided by a set of controllable generators, loads, or storage equipment potentially coupled with renewable energy sources. In this work, the DR flexibility is exclusively described in a power and energy domain, and it is modeled similarly to generic storage with a unitary efficiency. Then, additional constraints should be introduced to compute the "virtual state of charge" SOC*DR <sup>t</sup>* and keep it between the bounds (typically 0 and 100%) during the studied interval (8).

Two operating modes are envisioned for the designed DR. At first, when setting similar initial and final state-of-charge values (i.e., SOC*DR* <sup>0</sup> and SOC*DR <sup>t</sup>*=*t*) (typically 50%), this DR is managed in an "energy-shifting" mode, with an energy conservation constraint. The "energy-shifting" mode ensures the conservation of the consumed energy, especially through the constraints on the final values for the state of charge. Another operating mode called "energy shedding" is investigated in this paper. In "energy shedding" mode, the DR can be considered as a curtailable load (or an aggregation of curtailable loads). Within the proposed problem formulation, "energy shedding" is simply modeled by setting the "state of charge" at the beginning of the interval (SOC*DR* <sup>0</sup> ) to 100% and the SOC*DR <sup>t</sup>*=*<sup>t</sup>* at the end to 0% and with *PDR <sup>t</sup>* > 0 (i.e., equivalent to a DR discharge).

$$\begin{cases} SOC\_t^{DR} = SOC\_0^{DR} - \sum\_{k=1}^t \frac{100 \times P\_k^{DR} \times \Delta t}{E^{DRR}} \quad \forall t \in T\\ SOC^{DR} \le SOC\_t^{DR} \le \frac{\kappa}{SOC^{DR}}\\\\ SOC\_{t=T}^{DR} = SOC\_0^{DR} - \sum\_{k=1}^T \frac{100 \times P\_k^{DR} \times \Delta t}{E^{DRr}} \end{cases} \tag{8}$$

Multiplying the state-of-charge constraints on both the left and right hand side by the rated capacity *EDRr* allows removing the nonlinearity (introduced by the division of operating variable by the design variable). Thus, it allows us to solve the integrated management and sizing problem [55].

#### *3.3. Simplified Formulation for Thermal Constraints*

A first set of preliminary tests is intended to solve the proposed problem over a single simulated day using heuristics and meta-heuristics embedded in the MATLAB Optimization Toolbox (e.g., sequential quadratic programming, genetic algorithm). However, those runs suffered from expensive computational times and non-systematic convergence due to the problem size, especially with 1-min resolution for the reference temporal set. Additional complexity is introduced with the computation of the oil and hot-spot temperatures and aging as nonlinear constraints. In order to ensure the convergence within reasonable computational times so that different use cases can be investigated, the problem is linearized by taking advantage of the convexity of the thermal model.

Thus, a conventional piece-wise linearization (PWL) is implemented for the functions introduced in (1). The optimal breakpoints that minimize the linearization error are defined according to the method introduced in [56] and applied for different numbers of PWL segments *c* ∈ *C*. Figure 6a displays the obtained results when considering a three-block PWL for the nonlinear functions f<sup>1</sup> and f2. This PWL allows us to compute the slope coefficients *ac* for each function that is further used in the mathematical formulation illustrated for the function f<sup>1</sup> in Figure 6b.

**Figure 6.** Piece-wise linearization (PWL) process: (**a**) functions fitting and breakpoints; (**b**) mathematical formulation for function f1.

This typical formulation relies on the introduction of additional variables *δP*f1 *<sup>t</sup>*,*<sup>c</sup>* (similarly for f2) that represent the contribution of each block *c* at every time step *t*. Those semi-definite positive variables are bounded by the breakpoints previously defined, and their summed values should correspond to the transformer power *Ptr <sup>t</sup>* at every time step (9). Then, the functions are approximated with the contribution in each block weighted by the corresponding slope coefficient according to (10). Those estimations, for both functions f <sup>1</sup> and f2, are integrated in the equations for the transformer thermal model (2) and (3), which actually become constraints for the proposed optimization problem—i.e., mathematical constraints rather than physical constraints in Mathematical Language Programming. Since the approximated functions are convex, their PWL formulation should ultimately be introduced in the objective function so that the proper sequence of the activated block (i.e., *δP*f1 *<sup>t</sup>*,*c*, *δP*f2 *<sup>t</sup>*,*<sup>c</sup>* > 0) is ensured for different operating points at every time step. Thus, Equation (11) displays the updated version of the objective. Specific attention is given to the weight *αpwl* that implicitly performs an arbitrage between the "physical" objective, which is the minimization of the DR flexibility, and the "mathematical" objective that

allows the PWL for the convex functions considered. Especially, a high value of *αpwl* would tend to minimize those functions over the time horizon with temperatures that would remain much lower than the limits and could lead to oversized flexibilities. Preliminary runs and sensitivity analysis for various cases allowed setting a value of *αpwl* = 10−<sup>6</sup> so that the resulting DR flexibility is the minimum one. This allows us to reach the worst cases at the temperature limit but never at lower temperatures.

$$\begin{cases} \sum\_{\boldsymbol{c} \in \mathcal{C}} \delta p\_{t,\boldsymbol{c}}^{\mathrm{f1}} = P\_t^{\mathrm{tr}}\\ 0 \le \delta p\_{t,\boldsymbol{c}}^{\mathrm{f1}} \le p\_{\boldsymbol{c}}^{\mathrm{f1}} - p\_{\boldsymbol{c}-1}^{\mathrm{f1}} \end{cases} \quad \forall \boldsymbol{t} \in T \tag{9}$$

$$\mathbf{f}^{\mathbf{f}}(P\_t^{tr}) \leftarrow \mathbf{f}^{\mathbf{f}}\left(p\_0^{\mathbf{f}1}\right) + \sum\_{c \in C} a\_c^{\mathbf{f}1} \times \delta p\_{t,c}^{\mathbf{f}1} \quad \forall t \in T \tag{10}$$

$$obj: \min\left(P^{DRr} + \frac{E^{DRr}}{1\text{h}}\right) + a^{pwl} \times \sum\_{t \in T} \left(\begin{array}{c} \text{f}^1\left(p\_0^{t1}\right) + \sum\_{c \in \mathbb{C}} a\_c^{t1} \times \delta p\_{t,c}^{t1} \\ + \text{f}^1\left(p\_0^{t2}\right) + \sum\_{c \in \mathbb{C}} a\_c^{t2} \times \delta p\_{t,c}^{t2} \end{array}\right) \tag{11}$$

The optimization problem is written using the YALMIP library in MATLAB [57] and solved with CPLEX 12.0 (IBM, Armonk, New York, NY, USA). CPLEX is preferred here, since the linear programming in MATLAB is not as scalable over 100,000 variables as in CPLEX (in our problem 1 day ≈50,000 variables). Moreover, we observed that for long intervals, MATLAB *linprog* stops prematurely because it exceeds its allocated memory. As the consequence, we saw that on small horizons (one single day of our problem), CPLEX converges at least 25 times faster than MATLAB and continues converging when a number of variables grows at long intervals. As a reminder, the number of variables in the optimization problem becomes an issue for long horizons/reserve margins (see their relation in Figure 4) due to the 1-min time resolution required by the IEC standard [42] (up to 3 million variables for the longest horizons). The reader can see Section 3.5 on how the hour time resolution of *PDR <sup>t</sup>* was used for reducing the problem complexity while keeping temperatures (state variables) in 1-min time resolution in accordance with IEC standard [42].

Validation runs were performed without flexibility, and no thermal limits were identified. Figure 7 displays the results when comparing the output temperatures (top oil and hot spot) with the PWL approximation compared to the reference model.

**Figure 7.** PWL performances for different numbers of blocks—(**a**) top oil temperature; (**b**) hot spot temperature.

The approximation error is computed through the normalized **root-mean-square error** (in %) for different numbers of PWL blocks. This error decreases significantly with more blocks considered. A value of *C* = 6 is ultimately chosen as it allows us to reach a deviation of less than 1% with the reference.
