**Research on Multi-Time Scale Optimization Strategy of Cold-Thermal-Electric Integrated Energy System Considering Feasible Interval of System Load Rate**

#### **Bin Ouyang 1,2,\*, Zhichang Yuan 2, Chao Lu 2, Lu Qu <sup>2</sup> and Dongdong Li <sup>1</sup>**


Received: 15 July 2019; Accepted: 19 August 2019; Published: 22 August 2019

**Abstract:** The integrated energy system coupling multi-type energy production terminal to realize multi-energy complementarity and energy ladder utilization is of great significance to alleviate the existing energy production crisis and reduce environmental pollution. In this paper, the topology of the cold-thermal-electricity integrated energy system is built, and the decoupling method is adopted to analyze the feasible interval of load rate under the strong coupling condition, so as to ensure the "source-load" power balance of the system. Establishing a multi-objective optimization function with the lowest system economic operation and pollution gas emission, considering the attribute differences and energy scheduling characteristics of different energy sources of cold, heat and electricity, and adopting different time scales to optimize the operation of the three energy sources of cold, heat and electricity, wherein the operation periods of electric energy, heat energy and cold energy are respectively 15 min, 30 min and 1 h; The multi-objective problem is solved by standard linear weighting method. Finally, the mixed integer nonlinear programming model is calculated by LINGO solver. In the numerical simulation, the hotel summer front load parameters of Zhangjiakou, China are selected for simulation and compared with a single time scale system. The simulation results show that the multi-time scale system reduces the economic operation cost by 15.6% and the pollution gas emission by 22.3% compared with the single time scale system, it also has a wider feasible range of load rate, flexible time allocation, and complementary energy selection.

**Keywords:** cold-thermal-electricity integrated energy system; multi-energy complementation; feasible interval of load rate; multi-objective; multi time scales; optimize operation

#### **1. Introduction**

As the driving force for social development and the indispensable whole of daily life, energy plays an increasingly important role in modern life. As for the problems that the traditional fossil energy is exhausted, the energy utilization efficiency is low and the environmental pollution is serious which people have to seek an efficient, energy-saving, environment-friendly and renewable energy production and utilization mode. At the same time, due to the differences and limitations in the development of different energy systems, the planning, design, operation and control of each energy system are often separated and lack coordination with each other [1]. As a result, the overall energy utilization rate of the system is low, the energy supply of the system is safe, and the self-healing ability is not strong [2]. Therefore, Integrated Energy Systems (IES) [3–6] that couple multiple energy production and consumption modes, absorb a large amount of renewable energy, realize multi-energy complementation, and improve the overall energy utilization rate should emerge.

In recent years, many scholars have paid close attention to the research and application of integrated energy system. In Europe, the University of Manchester in the United Kingdom took the lead in developing a local integrated energy system, which integrates electricity/gas/thermal energy systems and a user interaction platform. The platform realizes three major functions: energy consumption mode, energy conservation strategy and demand-side response [2]; Denmark is the first country in the world to set a goal of completely separating from fossil fuels. It is expected to use 100% renewable energy by 2050, with emphasis on integrating different energy systems [7]; The University of Aachen and the German Federal Ministry of Economy and Environment launched the E-Energy Project [8] through demand-side response, achieving a high degree of integration of energy, information and capital, promoting the automation of energy service management, and successfully landing in Langenfel, Germany. The EU has set a target of carbon pollution from electricity production by 2050 [9], planned a new route for the EU power grid plan, and is committed to building a trans-European high-efficiency energy system by integrating the energy systems of various countries. In 2001, the U.S. Department of Energy put forward a comprehensive energy system development plan [10] aimed at improving the supply and utilization of clean and renewable energy and further improving the economy and reliability of the energy supply system. Japan established the Japan-Wireless City Alliance in 2010 [11] and is committed to the research of wireless city technology and the demonstration of the national integrated energy system. In China, a number of integrated energy system demonstration projects have been launched. Guangzhou Mingzhu Industrial Park, in combination with the future development direction and technical requirements of the city's power grid [12] that actively building an intelligent industrial park demonstration park for large-scale local consumption of renewable energy through cold/heat/electricity/gas multi-system coupling. Zhangbei wind/Light/Heat/Storage/Transmission Multi-energy Complementary Demonstration Project [13] in Zhangjiakou, Hebei province sets a precedent for large-scale multi-energy complementary power generation by comprehensively utilizing various energy storage and photo-thermal power generation technologies in order to build a strong smart grid. At the same time, the integrated energy system is also facing many difficulties and challenges in planning, modeling and optimizing operation.

In terms of optimal operation, document [14] establishes an economic optimal operation model of distributed energy system compatible with demand-side regulation, fully considers the cold/heat/electricity load in the distributed energy system which proposes quantum fireworks algorithm to solve the model. Literature [15] takes the cold/heat demand of 433 buildings as a model, compares the adaptability of two operation modes of heat determination and electricity determination in buildings, and finally points out that heat determination is not economical due to the high investment cost in the early stage, and the optimization of primary energy saving rate can improve the economic benefits of the factory. However, with the method of heat determination by electricity, when the demand for heat energy is strong, auxiliary equipment is needed to ensure the supply of heat energy. Excess heat is easily dissipated in the environment and wasted when the heat supply is excessive, so it is not ideal. In view of the irrationality of the methods of determining electricity by heat and determining heat by electricity, document [16] proposes energy scheduling algorithms aiming at minimizing operating cost, primary energy consumption (PEC) and carbon dioxide emission (CDE), respectively, points out that the three optimization algorithms do not have a common trend under common circumstances, and points out that the installation and operation of the system should be considered as long as PEC and CDE meet the standard requirements. However, using this method, different buildings need different analysis and evaluation to determine the optimal operating conditions. Literature [17] takes the minimum total operating cost as the optimization objective to maximize the utilization rate of resources. But, the problems mentioned above are all aimed at the single time scale operation of the system, and there are few researches on the multi-time scale characteristics of the integrated energy system. Documents [18] and [19] aim at the defects existing in the operation of the integrated energy system before the day, and use the two time scales of day-ahead load state estimation and real-time rolling correction to realize the optimization of real-time scheduling. Literature [20] considers the dynamic response characteristics of multi-time scale and multi-energy coupling systems, and starts with the modeling direction, constructs a network model for the combined calculation of electric/gas/heat/cold

multi-energy flows from both steady and dynamic aspects. In the process of solving the dynamic model, a hybrid stepping time domain simulation method for electromechanical transient simulation of power system and a medium and long-term transient simulation for non-power system are proposed. Finally, the electric heating coupling network is constructed for simulation verification. Document [21] proposes a multi-time scale multi-objective optimal joint scheduling model, which takes into account the different time scale characteristics of wind/water resources, electric/thermal load characteristics, peak load regulation capability, thermoelectric coupling characteristics of different thermal units, transmission capacity and system rotation standby requirements, and proposes a multi-time scale scheduling method and wind power prediction technology. In the short-term prediction technology, a day-to-day plan has been established to optimize the unit output of thermal generators, minimize operating costs and maximize the capacity to accommodate wind energy. Establish a time-of-day plan to optimize the power generation plans of all thermal, hydro and wind power plants on the basis of ultra-short-term wind power forecasting technology. However, in these studies, the attribute differences and energy dispatching characteristics of different energy sources are ignored, for example, the regulation of electric energy is more sensitive and its control mode is more flexible [22,23]; Thermal energy shows natural flexibility in regulation, with long response time and slow dynamic process [24]; However, natural gas needs to be added with conversion links during its use. On the premise that natural gas can be used in large quantities, it is inevitable to add large-capacity energy storage devices [25]; the cold energy change process is slow and requires a long dynamic time, which is not conducive to long-term storage [26]. At the same time, due to the strong coupling relationship of integrated energy systems, various energy systems influence each other, especially the non-linear devices of multiple coupled systems force the system to change frequently in operating load rate under the demand response, resulting in system "source-load" power imbalance. Therefore, based on the feasible interval conditions of system load rate, this paper studies the multi-time scale optimal operation of the cold-heat-electricity integrated energy system with the economic operation and the lowest pollutant gas emissions as multi-objective functions.

In the following content arrangement, Chapter 2 describes the topology and working principle of the cold-thermal-electricity integrated energy system. Chapter 3 analyzes the feasible range of load rate of integrated energy system. In Chapter 4, a multi-time scale optimization model is established with economic operation and minimum pollutant gas emissions as objective functions. In Chapter 5, the cold/heat/electricity front load of Jianguo hotel in Zhangjiakou, China is selected for simulation analysis and compared with a single time scale system. Chapter 6 summarizes this article.

#### **2. Topological Structure and Working Principle of Cold-Thermal-Electricity Integrated Energy System**

The integrated energy system has various structures, different forms and complicated coupling relationships, there are various energy production and conversion devices in the system. In this paper, the topology structure adopted for the research on the optimal operation of the cold-thermal-electricity integrated energy system is shown in Figure 1.

The cold-thermal-electricity integrated energy system is micro-energy network level. The equipment in the system mainly includes: gas internal combustion engine (GE), flue gas absorption heat pump (AP), cylinder liner water heat exchanger (JW), absorption refrigerator (AC), electric boiler (EB) and electric refrigerator (EC), and two kinds of energy storage equipment including electricity storage (ST) and heat storage (HS) are added. At the same time, in order to make full use of local reliable solar energy resources, a photovoltaic generator set (PV) is added and connected to the power network to ensure the balance of power supply and demand in the system. The whole system takes a gas internal combustion engine and a flue gas absorption heat pump as the core, the gas internal combustion engine directly supplies part of the electric load by consuming natural gas and generating electric energy, and the high-temperature steam generated during working is converted into hot water to supply the heat load through a cylinder liner water heat exchanger; The flue gas generated during

the combustion of natural gas can be absorbed and utilized by most of the flue gas absorption heat pumps and converted into heat energy and cold energy for direct supply to users; The absorption refrigerator converts part of the heat energy on the heat bus into cold energy to supply the cold load for use; When the system needs more heat energy or cold energy, some of the heat energy and cold energy deficiency can be made up through the work of the electric boiler or electric refrigerator. The system is also connected with two energy storage devices of heat storage and electricity storage to ensure that the system has sufficient margin of electric/thermal power capacity and increase the stability of the system. In addition, the active connection of photovoltaic generator sets not only effectively utilizes solar energy resources, but also increases the environmental protection and economic benefits of the system. When the demand for electric energy is large, the system can interact with the power grid. At the same time, in order to reduce the construction cost and coordination cost of the information channel and physical channel between the system and the power grid, the system adopts the principle of "connecting to the grid without power output " to purchase electric energy from the power grid, so as to make up for the shortage of electric energy in the system and ensure the stable operation of the system.

**Figure 1.** Topology diagram of cold-thermal-electricity integrated energy system.

#### **3. The Feasible Range of System Load Rate**

The integrated energy system realizes the coupling relation of multiple types of heterogeneous energy flows, and the system has strong coupling, nonlinearity, multi-energy complementation and mutual influence characteristics. For the cold-thermal-electricity integrated energy system, even if the operating characteristics and operating condition change parameters of each equipment in the system are clearly defined, there are many uncertainties in the energy output of the system under the condition that various energy input terminals are determined. Especially in the research process, it is found that when the cold-thermal-electricity integrated energy system is directly connected to the load side, the output power of the system does not match the load and the load rate of the system cannot find a reliable operation interval, so it cannot operate stably. Therefore, it is necessary to analyze the feasible range of system load rate.

Because of the close coupling relationship between the cold-thermal-electricity integrated energy system, it is difficult to directly analyze the load rate of the system. Therefore, it is possible to decouple the system and analyze the operational load rate characteristics of each decoupling subsystem, thus ensuring the "source-load" power balance of the system.

The decoupling method in this paper is as follows:


Therefore, the feasible range of the load rate of each decoupling subsystem of the cold-thermalelectricity integrated energy system is shown in the following Figure 2:

**Figure 2.** Load rate feasible interval of decoupling subsystem. (**a**) Decoupling power subsystem; (**b**) Decoupling the low load characteristics of the thermal energy subsystem; (**c**) Decoupling the high load characteristics of the thermal energy subsystem; (**d**) Decoupling the cold energy subsystem.

Since the loads in the cold-thermal-electricity integrated energy system are independent of each other and do not interfere with each other. For the analysis of the feasible range of system load rate, as shown in Figure 2a, when the electric power load is *P*0, it intersects with the upper and lower limits of electric power at two points of load rate a*Pi* and a*Pj*, then the feasible load rate of the power subsystem will be between (a*Pi*, a*Pj*), and the electric power output power is shown by shadow *S*<sup>1</sup> in the figure. There is a critical jump *QH* in the upper limit of the power output of the decoupled thermal energy subsystem, where the system load rate is a*H*. When the thermal energy load received by the system is lower than *QH*, it is called the low load characteristic of the thermal energy subsystem, and when the thermal energy load received by the system is higher than *QH*, it is called the high load characteristic of the thermal energy subsystem. As shown in Figure 2b, under the condition of low load characteristics, the upper and lower limits of the thermal energy load *QH*<sup>0</sup> of the system and the power output of the thermal energy decoupling subsystem intersect at two points of the load rates a*Hi* and a*Hj*, so that the system load rate is limited between (a*Hi*, a*Hj*), and the thermal energy output range of the thermal energy subsystem is shaded *S*<sup>2</sup> in Figure 2b; When the thermal energy load *QH*<sup>1</sup> is higher than *QH*, the thermal energy subsystem exhibits a high load characteristic, as shown in Figure 2c, the load rate of the system will be higher than a*Hk*, and the thermal energy output range is shown in shadow *S*3. For the cold energy subsystem, as shown in Figure 2d, the lower limit of the cold energy output power is stable, while the upper limit of the output power varies approximately linearly, so that the cold energy system has a wider load rate selection range when outputting power outward. When the system cold energy load is *QC*0, the system load rate can fall between (a*Ci*, 1), and the cold energy output is shown in shadow *S*4. To sum up, the analysis of the operating characteristics of the load rate of each decoupling subsystem shows that the loads of inter-cold, heat and electricity in the system are independent and uncertain. In order to balance the "source-load" power of the system, the system load rate A should be in the same interval as the operational load rate of each decoupling subsystem, namely:

$$\begin{array}{l} a \in (a\_{\text{Pi}\prime}, a\_{Pj}) \cap (a\_{\text{Hi}\prime}, a\_{Hj}) \cap (a\_{\text{Ci}\prime}\mathbf{1})\\\text{or} \ a \in (a\_{\text{Pi}\prime}, a\_{Pj}) \cap (a\_{\text{Hk}\prime}\mathbf{1}) \cap (a\_{\text{Ci}\prime}\mathbf{1}) \end{array} \tag{1}$$

However, when the load rate intervals of each decoupling subsystem have no intersection, the "source-load" power of the system is unbalanced and cannot operate stably. Considering that the system load fluctuates within a stable range, it is possible to change the operating conditions of the system by adjusting the system equipment parameters and widen the intersection interval of load rates so that the system has a stable feasible interval. This paper will also study the optimal operation of the cold-thermal-electricity integrated energy system on the basis of the feasible load rate range of the system.

#### **4. Optimize Operation**

The operation process of the cold-thermal-electricity integrated energy system is complex, with many parameters and various energy outputs. This paper optimizes the operation of the cold-thermal-electricity integrated energy system with multi-objective and multi-time scales. In terms of multi-objective function, the system has the lowest operating cost and the least pollutant gas emission. In terms of multi-time scale, considering the difference of attributes of different energy sources and the different flexible characteristics of operation, for example, short power dispatching period, flexible regulation and fast power release and absorption, therefore, power is operated at Δ*t*2,i = 15 min as one dispatching time; However, heat energy is easier to adjust and store than cold energy. Therefore, heat energy is operated at Δ*t*2,j = 30 min as a scheduling time, while cold energy is operated at Δ*t*<sup>1</sup> = 1 h as a scheduling time. In this way, the operation periods of cold, heat and electricity are multiple of each other, which is beneficial for multiple energy sources to supplement each other and ensures the consistency of multi-time scale coordination. At the same time, in order to reduce system losses and adjustment costs, the load rate of the system is kept unchanged for one hour, i.e., 4 periods of power adjustment.

#### *4.1. Objective Function*

The objective function is to construct a multi-objective function with the lowest total operating cost of the system and the lowest pollutant gas emissions, as follows:

$$\min \left\{ \begin{array}{l} F\_{\text{run}} = \sum\_{t\_1=1}^{24} \sum\_{i=1}^{4} \left[ F\_{\text{grid}}(t\_1, t\_2, i) + F\_{\text{gas}}(t\_1, t\_2, i) + F\_{\text{min}}(t\_1) \right] \\\ F\_{\text{pull}} = \sum\_{t\_1=1}^{24} \sum\_{i=1}^{4} \left( a\_{\text{surr}} \cdot P\_{\text{grid}}(t\_1, t\_2, i) + a\_{\text{trun}} \cdot P\_{\text{grid}}(t\_1, t\_2, i) + a\_{\text{PGE}} \cdot P\_{\text{CE}}(t\_1, t\_2, i) \right) \end{array} \tag{2}$$

In the total operating cost, the electricity purchase cost, natural gas purchase cost and system equipment maintenance cost of the system and power grid are considered; Pollution gas emissions from the power supply side of the power grid, power loss from the power transmission lines of the power grid, and pollution gas emissions from the operation of the gas internal combustion engine are considered in the pollution gas emissions.

Where *Frun* is the total operating cost of the system; *Fgrid*(*t*1,*t*2,i) is the electricity purchase cost of the system and the power grid; *Fgas*(*t*1,*t*2,i) purchase natural gas for the system; *Fmain*(*t*1) is the maintenance cost of system equipment; *F*poll is the amount of pollutant gas discharged; α*sour* is the pollutant gas emission coefficient on the power supply side of the power grid; α*tran* is the pollutant gas emission coefficient delivered by power grid lines; *Pgrid*(*t*1,*t*2,i) purchases power for the system and power grid; α*PGE* is the pollutant emission coefficient of gas internal combustion engine. *PGE*(*t*1,*t*2,i) is the output electric power of the gas internal combustion engine; *t*<sup>1</sup> is expressed as the number of hours in a day. *t*<sup>2</sup> is expressed as the number of minutes in an hour. Since electric energy runs at Δ*t*2,i = 15 min as a scheduling time, electric energy is scheduled 4 times in an hour, *i* = 1, 2, 3 and 4. Similarly, heat energy is dispatched twice in one hour, with *j* = 1 and 2.

Among them, the system electricity purchase cost is specifically expressed as follows:

$$F\_{\rm grid}(t\_1, t\_{2,i}) = P\_{\rm grid}(t\_1, t\_{2,i}) \cdot \Delta t\_{2,i} \cdot f\_{\rm grid}(t\_1, t\_{2,i}) \tag{3}$$

where, *fgrid*(*t*1,*t*2,i) is the real-time electricity price of the power grid.

The system purchase cost of natural gas is specifically expressed as follows:

$$F\_{\rm gus}(t\_1, t\_{2,i}) = V\_{\rm gus}(t\_1, t\_{2,i}) \cdot \Delta t\_{2,i} \cdot f\_{\rm gus}(t\_1, t\_{2,i}) \tag{4}$$

where *Vgas*(*t*1,*t*2,i) is the volume of natural gas consumed by the system; *fgas*(*t*1,*t*2,i) is the price of natural gas.

The system equipment maintenance costs are specified as follows:

*Fmain*(*t*1) <sup>=</sup> <sup>4</sup> *i*=1 2 *j*=1 ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ *kGE*[*PGE*(*t*1, *t*2,*i*)] · Δ*t*2,*<sup>i</sup>* · *PGE*(*t*1, *t*2,*i*) + *kPV*[*PPV*(*t*1, *t*2,*i*)] · Δ*t*2,*<sup>i</sup>* · *PPV*(*t*1, *t*2,*i*) + *kAP*.*cool*[*QAP*.*cool*(*t*1, *t*2,*i*)] · Δ*t*2,*<sup>i</sup>* · *QAP*.*cool*(*t*1,*t*2,*i*)+ *kAP*.*heat*[*QAP*.*heat*(*t*1, *t*2,*i*)] · Δ*t*2,*<sup>i</sup>* · *QAP*.*heat*(*t*1, *t*2,*i*) + *kAC*.*heat*[*QAC*.*heat*(*t*1, *t*2,*<sup>j</sup>* )] · Δ*t*2,*<sup>j</sup>* · *QAC*.*heat*(*t*1, *t*2,*<sup>j</sup>* )+ *kbatt*.*dis*/*cha*[*Pbatt*.*dis*/*cha*(*t*1, *<sup>t</sup>*2,*i*)] · <sup>Δ</sup>*t*2,*<sup>i</sup>* · *Pbatt*.*dis*/*cha*(*t*1, *<sup>t</sup>*2,*i*) + *kstor*.*dis*/*cha*[*Qstor*.*dis*/*cha*(*t*1, *<sup>t</sup>*2,*<sup>j</sup>* )] · <sup>Δ</sup>*t*2,*<sup>j</sup>* · *Qstor*.*dis*/*cha*(*t*1,*t*2,*<sup>j</sup>* ) ⎫ ⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎭ (5)

where *kGE*[*PGE*(*t*1,*t*2,i)] is the maintenance coefficient of the gas internal combustion engine at different output powers; *kPV*[*PPV*(*t*1,*t*2,i)] is the maintenance coefficient of photovoltaic generator set; *PPV*(*t*1,*t*2,i) is the generating power of photovoltaic generator set; *kAP.cool*[*QAP.cool*(*t*1,*t*2,i)] is the cold power maintenance coefficient of the flue gas absorption heat pump; *QAP.heat*(*t*1,*t*2,i) is the cold power output by the flue gas absorption heat pump; *kAC.heat*[*QAC.heat*(*t*1,*t*2,j)] is the thermal power maintenance coefficient of the flue gas absorption heat pump; *QAP.heat*(*t*1,*t*2,i) is the output heat power of the flue gas absorption heat pump; *kAC.heat*[*QAC.heat*(*t*1,*t*2,j)] is the maintenance coefficient of absorption chiller; *QAC.heat*(*t*1,*t*2,j) is the heat power absorbed by the absorption refrigerator; *kbatt.dis*/*cha*[*Pbatt.dis*/*cha*(*t*1,*t*2,i)], *kstor.dis*/*cha*[*Qstor.dis*/*cha*(*t*1,*t*2,j)] are the power maintenance coefficients of power storage equipment and heat storage equipment respectively; *Pbatt.dis*/*cha*(*t*1,*t*2,i) and *Qstor.dis*/*cha*(*t*1,*t*2,j) are the interactive power of electric storage equipment and heat storage equipment.

#### *4.2. Constraints*

The cold-thermal-electric integrated energy system constraints include: equipment model constraints and power balance constraints.

#### 4.2.1. Equipment Model Constraints

Among them, the equipment model includes a gas internal combustion engine model [27,28], a flue gas absorption heat pump model [29–31], an absorption refrigerator model [32], a cylinder liner water heat exchanger model, an electric boiler model, an electric refrigerator model, a photovoltaic generator set model [33], a heat and electricity storage model [34–36].

#### • Gas engine model.

Among them, gas-fired internal combustion engines have good electrical energy and thermal energy output characteristics which can be divided into power generation, heat generation and consumption of natural gas.

Power generation of gas internal combustion engine, the maximum real-time output power of an internal combustion engine is limited to *P*max:

$$\begin{cases} \eta\_{\rm GE.clx}(t\_1, t\_2; 1) = a\_3 \left( \frac{P\_{\rm GE}(t\_1, t\_{2,i})}{P\_{\rm max}} \right)^3 + a\_2 \left( \frac{P\_{\rm GE}(t\_1, t\_{2,i})}{P\_{\rm max}} \right)^2 + a\_1 \left( \frac{P\_{\rm GE}(t\_1, t\_{2,i})}{P\_{\rm max}} \right) + a\_0\\ \qquad P\_{\rm GE}(t\_1, t\_{2,1}) = P\_{\rm GE}(t\_1, t\_{2,2}) = P\_{\rm GE}(t\_1, t\_{2,3}) = P\_{\rm GE}(t\_1, t\_{2,4})\\ \qquad \left| P\_{\rm GE}(t\_1, t\_{2,i}) - P\_{\rm GE}(t\_1, t\_{2,i} - 1) \right| \le P\_{\rm GE.\max} \end{cases} \tag{6}$$

Heating power of gas internal combustion engine:

$$\begin{cases} \begin{array}{c} Q\_{\rm GE.heat}(t\_1, t\_{2,i}) = \frac{P\_{\rm GE}(t\_1, t\_{2,i})}{\eta\_{\rm GE.dev}(t\_1, t\_{2,i})} (1 - \eta\_{\rm GE.clec}(t\_1, t\_{2,i}) - \eta\_{\rm L})\\\ Q\_{\rm GE.heat}(t\_1, t\_{2,1}) = Q\_{\rm GE.heat}(t\_1, t\_{2,2}) = Q\_{\rm GE.heat}(t\_1, t\_{2,3}) = Q\_{\rm GE.heat}(t\_1, t\_{2,4}) \end{array} \end{cases} \tag{7}$$

Natural gas consumed by gas internal combustion engines:

$$V\_{\rm gas}(t\_1, t\_{2,i}) = \frac{P\_{\rm GE}(t\_1, t\_{2,i}) \cdot \Delta t\_2}{\eta\_{\rm gas} \cdot \eta\_{\rm GE, clec}(t\_1, t\_{2,i}) \cdot LHV} \tag{8}$$

In the formula, the output electric power *PGE*(*t*1,*t*2,i) of the gas internal combustion engine remains unchanged for one hour, i.e., four periods of electric energy output; η*GE.elec*(*t*1,*t*2,i) is the power generation efficiency of the gas internal combustion engine; *P*max is the rated power of gas internal combustion engine. The thermal power *QGE.heat*(*t*1,*t*2,i) output by the gas internal combustion engine remains constant within one hour. η*<sup>L</sup>* is the inherent loss rate of gas internal combustion engine; *PGE.*max is the output gradient constraint of gas internal combustion engine. LHV is the low calorific value of natural gas; η*gas* is the natural gas utilization rate of gas internal combustion engines; *a*3, *a*2, *a*<sup>1</sup> and *a*<sup>0</sup> are fitting constants respectively.

• Smoke absorption heat pump model.

Smoke absorption heat pump model can be divided into three parts: absorbing smoke, outputting heat energy and cold energy.

Smoke absorption heat pump absorbs smoke, absorption heat pump absorbed the flue gas temperature are directly affected parameter *COPAP*:

$$\begin{cases} \begin{aligned} T(t\_1, t\_{2,i}) &= b\_1 \cdot (\frac{P\_{GE}(t\_1, t\_{2,i})}{P\_{\max}}) + b\_0 \\ \mathcal{C}\_w(t\_1, t\_{2,i}) &= b\_3 \cdot T(t\_1, t\_{2,i}) + b\_2 \\ \mathcal{C}OP\_{AP}(t\_1, t\_{2,i}) &= b\_5 \cdot (\frac{P\_{GE}(t\_1, t\_{2,i})}{P\_{\max}}) \\ \lambda\_{\text{hat}}(t\_1, t\_{2,i}) + \lambda\_{\text{cool}}(t\_1, t\_{2,i}) &= 1 \end{aligned} \tag{9}$$

⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩

> ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩

At the same time, the real-time power of flue gas absorption heat pump has linear constraint values *QAP.heat.*max and *QAP.cool.*max.

Smoke absorbs heat output from heat pump:

$$\begin{aligned} Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,l}) &= \lambda\_{\text{hart}}(t\_1, t\_{2,l}) \cdot \mathbb{C}\_{w}(t\_1, t\_{2,l}) \cdot (T(t\_1, t\_{2,l}) - T\_{\text{hart}}) \cdot \text{COP}\_{AF}(t\_1, t\_{2,l}) \cdot L\_{\text{hart}}(t\_1, t\_{2,l}) \cdot \eta\_{AP\text{-}h\text{-}t} \\ Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,1}) &= Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,2}) = Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,3}) = Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,4}) \\ & \left| Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,l}) - Q\_{AP\text{-}h\text{-}t}(t\_1, t\_{2,l} - 1) \right| \le Q\_{AP\text{-}h\text{-}t, \text{max}} \\ & 0 \le L\_{\text{hart}}(t\_1, t\_{2,l}) \le L\_{\text{hart}, \text{max}} \end{aligned} \tag{10}$$

Smoke absorption heat pump outputs cold energy:

$$\begin{array}{c} Q\_{AP\,\text{coll}}(t\_1, t\_{2,l}) = \lambda\_{\text{coll}}(t\_1, t\_{2,l}) \cdot \mathbb{C}\_{\text{-}0}(t\_1, t\_{2,l}) \cdot (T(t\_1, t\_{2,l}) - T\_{\text{coll}}) \cdot \mathbb{C}OP\_{AP}(t\_1, t\_{2,l}) \cdot L\_{\text{coll}}(t\_1, t\_{2,l}) \cdot \eta\_{AP\,\text{coll}} \\ \mathbb{Q}\_{AP\,\text{coll}}(t\_1, t\_{2,l}) = \mathbb{Q}\_{AP\,\text{coll}}(t\_1, t\_{2,l}) = \mathbb{Q}\_{AP\,\text{coll}}(t\_1, t\_{2,l}) = \mathbb{Q}\_{AP\,\text{coll}}(t\_1, t\_{2,l}) \\ \left[ Q\_{AP\,\text{coll}}(t\_1, t\_{2,l}) - Q\_{AP\,\text{coll}}(t\_1, t\_{2,l} - 1) \right] \le Q\_{AP\,\text{coll}\,\text{max}} \\ 0 \le L\_{\text{val}}(t\_1, t\_{2,l}) \le L\_{\text{val}} \,\text{max} \end{array} \tag{11}$$

where in *T*(*t*1,*t*2,i) is the inlet temperature of the flue gas absorption heat pump; *CW*(*t*1,*t*2,i) is the specific heat capacity of hot water at different temperatures; *COPAP*(*t*1,*t*2,i) is the energy efficiency coefficient of flue gas absorption heat pump; The heating power *QAP.heat*(*t*1,*t*2,i)) and the cooling power *QAP.cool*(*t*1,*t*2,i) of the flue gas absorption heat pump remain unchanged within one hour. λ*heat*(*t*1,*t*2,i) and λ*cool*(*t*1,*t*2,i) are the heating ratio and cooling ratio of flue gas of the flue gas absorption heat pump respectively; *Theat* and *Tcool* are hot water outlet temperature and cold water outlet temperature respectively. *Lheat*(*t*1,*t*2,i) and *Lcool*(*t*1,*t*2,i) are the hot water and cold water flows of the flue gas absorption heat pump respectively; *Lheat.*max and *Lcool.*max are the maximum heating and cooling flows respectively; ηAP.heat and ηAP.cool are the heating and cooling efficiency of flue gas absorption heat pump respectively. *QAP.heat.*max is the gradient constraint of heating power output of flue gas absorption heat pump; *QAP.cool.*max is the gradient constraint of cooling power output of flue gas absorption heat pump; *b*5, *b*4, *b*3, *b*2, *b*<sup>1</sup> and *b*<sup>0</sup> are fitting constants respectively.

• Absorption refrigerator model.

$$\begin{cases} Q\_{AC\,\text{cool}}(t\_1) = \text{COP}\_{AC} \cdot \sum\_{j=1}^{2} Q\_{AC\,\text{heat}}(t\_1, t\_{2,j})\\ Q\_{AC\,\text{heat},\text{min}} \le Q\_{AC\,\text{heat}}(t\_1, t\_{2,j}) \le Q\_{AC\,\text{heat},\text{max}}\\ \left| Q\_{AC\,\text{cool}}(t\_1) - Q\_{AC\,\text{cool}}(t\_1 - 1) \right| \le Q\_{AC\,\text{cool},\text{max}} \end{cases} \tag{12}$$

where *QAC.cool*(*t*1) is the cold power output by the absorption refrigerator; *COPAC* is the energy efficiency coefficient of absorption refrigerator; *QAC.heat.*min and *QAC.heat.*max are the minimum and maximum thermal power absorbed by the absorption chiller respectively. *QAC.cool.*max is the output gradient constraint of absorption chiller.

In order to stabilize the output power of absorption refrigerator, the absolute value of the difference between the next output power and the previous one is *QAC.cool.*max.

• Cylinder liner water heat exchanger model.

$$Q\_{\rm JW}(t\_1, t\_{2,j}) = \eta\_{\rm JW} \cdot \sum\_{i=2j-1}^{2j} Q\_{\rm GE}(t\_1, t\_{2,i}) \tag{13}$$

where *QJW*(*t*1,*t*2,*j*) is the output thermal power of the cylinder liner water heat exchanger; η*JW* is the heat transfer efficiency of cylinder liner water heat exchanger.

• Electric boiler model.

$$\begin{cases} Q\_{EB}(t\_1, t\_{2,j}) = COP\_{EB} \cdot \sum\_{i=2j-1}^{2j} P\_{EB}(t\_1, t\_{2,i}) \\ P\_{EB, \text{min}} \le P\_{EB}(t\_1, t\_{2,i}) \le P\_{EB, \text{max}} \\ \left| Q\_{EB}(t\_1, t\_{2,j}) - Q\_{EB}(t\_1, t\_{2,j} - 1) \right| \le Q\_{EB, \text{max}} \end{cases} \tag{14}$$

In the formula, *PEB*(*t*1,*t*2,i) is the input electric power of the electric boiler; *QEB*(*t*1,*t*2,j) is the output thermal power of the electric boiler; *COPEB* is the energy production coefficient of electric boilers; *PEB*.min and *PEB*.max are respectively the minimum and maximum electric power of electric boilers; *QEB*.max is the output gradient constraint of the electric boiler.

• Electric refrigerator model.

$$\begin{cases} Q\_{EC}(t\_1) = COP\_{EC} \cdot \sum\_{i=1}^{4} P\_{EC}(t\_1, t\_{2,i}) \\ P\_{EC, \text{min}} \le P\_{EC}(t\_1, t\_{2,i}) \le P\_{EC, \text{max}} \\ \left| Q\_{EC}(t\_1) - Q\_{EC}(t\_1 - 1) \right| \le Q\_{EC, \text{max}} \end{cases} \tag{15}$$

where *PEC*(*t*1,*t*2,i) is the input electric power of the electric refrigerator; *QEC*(*t*1) is the output cold power of the electric refrigerator; *COPEC* is the energy efficiency coefficient of electric refrigerator. *PEC*.min and *PEC*.max are the minimum and maximum electric power of electric refrigerator respectively. *QEC*.max is the output gradient constraint of the electric refrigerator.

• Photovoltaic generator set model.

$$P\_{PV}(t\_1, t\_{2,i}) = P\_{STC} \frac{G\_{INC}(t\_1, t\_{2,i})}{G\_{STC}} [1 - k(T\_{out}(t\_1, t\_{2,i}) - T\_s)] \tag{16}$$

where, *PSTC* is the rated output of photovoltaic generator set; *GING*(*t*1,*t*2,i) is the real-time irradiation intensity; *GSTC* is the rated irradiation intensity of photovoltaic generator set; *K* is the power generation coefficient of photovoltaic generator set; *Tout*(*t*1,*t*2,i) is the external temperature; *T*<sup>s</sup> is the reference temperature of the generator set.

• Electricity storage model.

⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩

In the model of electricity storage model, the power balance and charge and discharge constraints are considered.

$$\begin{aligned} E\_{\text{batt}}(t\_1, t\_{2,i}) &= (1 - k\_L) \cdot E\_{\text{batt}}(t\_1, t\_{2,i} - 1) + \left[ \eta\_{\text{batt.cdu}} \cdot P\_{\text{batt.cdu}}(t\_1, t\_{2,i}) - \frac{P\_{\text{batt.cdu}}(t\_1, t\_{2,i})}{\eta\_{\text{batt.cdu}}} \right] \cdot \Delta t\_{2,i} \\ P\_{\text{batt.cdu.min}} &\le P\_{\text{batt.cdu}}(t\_1, t\_{2,i}) \le P\_{\text{batt.cdu.max}} \\ P\_{\text{batt.dds.min}} &\le P\_{\text{batt.dds}}(t\_1, t\_{2,i}) \le P\_{\text{batt.dds.cux}} \\ E\_{\text{batt.m.m}} &\le E\_{\text{batt}}(t\_1, t\_{2,i}) \le E\_{\text{batt.cux}} \end{aligned} \tag{17}$$

where in *Ebatt*(*t*1,*t*2,i) is that real-time capacity of the pow storage equipment; *kL* is the self-loss coefficient of power storage equipment; η*batt.cha* and η*batt.dis* are the charging and discharging efficiency of power storage equipment respectively. *Pbatt.cha*(*t*1,*t*2,i) and *Pbatt.dis*(*t*1,*t*2,i) are the charging and discharging power of the power storage equipment respectively. *Pbatt.dis.*max and *Pbatt.dis.*min are the maximum and minimum discharge powers of power storage equipment respectively. *Pbatt.cha.*max and *Pbatt.cha.*min are respectively the maximum and minimum charging power of power storage equipment. *Ebatt*.max and *Ebatt*.min are the maximum and minimum storage capacities of power storage equipment respectively. ⎧ ⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎩

#### • Thermal storage model.

Thermal storage model have similar model constraints to electricity storage model.

$$\begin{aligned} B\_{\text{star}}(t\_1, t\_{2,j}) &= (1 - k\_s) \cdot B\_{\text{star}}(t\_1, t\_{2,j} - 1) + [\eta\_{\text{star}.\text{clar}} \cdot Q\_{\text{star}.\text{clar}.\text{clar}}(t\_1, t\_{2,j}) - \frac{Q\_{\text{star}.\text{clar}.\text{clar}}(t\_1, t\_{2,j})}{\eta\_{\text{star}.\text{clar}.\text{}}}] \cdot \Delta t\_{2,j} \\ Q\_{\text{star}.\text{clar}.\text{min}} &\leq Q\_{\text{star}.\text{clar}.\text{clar}}(t\_1, t\_{2,j}) \leq Q\_{\text{star}.\text{clar}.\text{max}} \\ Q\_{\text{star}.\text{clar}.\text{min}} &\leq Q\_{\text{star}.\text{clar}.\text{}}(t\_1, t\_{2,j}) \leq Q\_{\text{star}.\text{clar}.\text{max}} \\ B\_{\text{star}.\text{min}} &\leq S\_{\text{star}}(t\_1, t\_{2,j}) \leq B\_{\text{star}.\text{max}} \end{aligned} \tag{18}$$

where *Bstor*(*t*1,*t*2,j) is the real-time capacity of heat storage equipment; *ks* is the self-loss coefficient of heat storage equipment; η*stor.cha* and η*stor.dis* are the heat absorption and heat release efficiency of heat storage equipment respectively. *Qstor.cha*(*t*1,*t*2,j) and *Qstor.dis*(*t*1,*t*2,j) are the heat absorption and heat release power of heat storage equipment respectively. *Qstor.cha*.max and *Qstor.cha*.min are the maximum and minimum heat absorption powers of heat storage equipment respectively. *Qstor.dis*.max and *Qstor.dis.*min are the maximum and minimum heat release power of heat storage equipment respectively. *Bstor*.max and *Bstor*.min are the maximum and minimum capacity constraints of heat storage equipment respectively.

#### 4.2.2. Power Balance Constraints

The power balance constraints of cold, heat and electricity are satisfied in the system.

• Electric power balance constraint.

$$\begin{aligned} P\_{\text{grid}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) + P\_{\text{PV}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) &\quad + P\_{\text{CE}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) + P\_{\text{hatt.dis}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) \cdot D\_{\text{hatt.dis}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) = \\ P\_{\text{hatt.cha}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) \cdot D\_{\text{hatt.cha}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) &+ P\_{\text{elc}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) + P\_{\text{EB}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) + P\_{\text{EC}}(\mathbf{t}\_{1},\mathbf{t}\_{2\perp}) \end{aligned} \tag{19}$$

In the formula, *Dbatt.dis*(*t*1,*t*2,i) and *Dbatt.cha*(*t*1,*t*2,i) are respectively discharge and charge variables of power storage equipment; *Pele*(*t*1,*t*2,i) is the power load.

• Thermal power balance constraint.

$$\begin{aligned} \underset{\text{(2)}}{Q\_{\text{PW}}(t\_1, t\_{2,j})} + \sum\_{i=2j-1}^{2j} Q\_{\text{AP-ball}}(t\_1, t\_{2,i}) &\quad + Q\_{\text{EE}}(t\_1, t\_{2,j}) + Q\_{\text{site-dis}}(t\_1, t\_{2,j}) \cdot D\_{\text{site-dis}}(t\_1, t\_{2,j}) = \\ &\quad Q\_{\text{site-dis}}(t\_1, t\_{2,j}) \cdot D\_{\text{site-dis}}(t\_1, t\_{2,j}) + Q\_{\text{heat}}(t\_1, t\_{2,j}) + Q\_{\text{Att-hant}}(t\_1, t\_{2,j}) \end{aligned} \tag{20}$$

In the formula, *Dstor.dis*(*t*1,*t*2,j) and *Dstor.cha*(*t*1,*t*2,j) are the heat release and heat absorption variables of the heat storage equipment respectively; *Qheat*(*t*1,*t*2,j) is thermal load.

• Cold power balance constraint.

$$Q\_{AC\text{-}col}(t\_1) + Q\_{EC}(t\_1) + \sum\_{i=1}^{4} \left[ Q\_{AP\text{-}col}(t\_1, t\_{2,i}) \right] = Q\_{coul}(t\_1) \tag{21}$$

where *Qcool*(*t*1)is the cooling load.

#### *4.3. Solution*

#### 4.3.1. Multi-Objective Solution Method

This model is a multi-objective mixed integer nonlinear programming model. The multi-objective problem is transformed into a single-objective problem that is easy to solve by adopting a scalar linear weighting method. The results under different weight conditions are compared by changing the weight coefficient to obtain the optimal result. The scalar process is as follows:

*Energies* **2019**, *12*, 3233

Firstly, the optimal values of economic operation and pollutant gas emission under single-objective conditions are solved respectively, and *Frun*.min and *Fpoll*.min are obtained.

$$\min \quad F\_{\text{run}} = \sum\_{t\_1=1}^{24} \sum\_{i=1}^{4} \left\{ F\_{\text{grid}}(t\_1, t\_{2,i}) + F\_{\text{gas}}(t\_1, t\_{2,i}) + F\_{\text{rain}}(t\_1) \right\} \tag{22}$$

$$\min \quad F\_{\text{poll}} = \sum\_{t\_1=1}^{24} \sum\_{i=1}^{4} \left\{ \alpha\_{\text{ssurr}} \cdot P\_{\text{grid}}(t\_1, t\_{2,i}) + \alpha\_{\text{tran}} \cdot P\_{\text{grid}}(t\_1, t\_{2,i}) + \alpha\_{\text{PGE}} \cdot P\_{\text{GE}}(t\_1, t\_{2,i}) \right\} \tag{23}$$

Secondly, the multi-objective optimization problem is transformed into a single-objective problem solving calculation through a linear weighting method. The solving process is as follows:

$$\min \, F = k\_{run} \cdot \frac{F\_{run}}{F\_{run.min}} + k\_{pol} \cdot \frac{F\_{poll}}{F\_{poll.min}}, \, k\_{run} + k\_{pull} = 1 \tag{24}$$

where *F* is the mixed objective function value; *krun* is the economic operation weight coefficient; *kpoll* is the weight coefficient of pollutant gas emission.

By changing the values of the weight coefficients *krun* and *kpoll* to adjustment and optimization results of the objective function are compared, the sum of weight coefficients is always 1. In this paper, 0.1 is selected as the discrimination degree of two parameters. Therefore, the value of the weight coefficient (*krun*, *kpoll*) = (0.1,0.9), (0.2,0.8) ... (0.9,0.1) changes.

Finally, by changing different weight coefficients, the optimal operating conditions under different weights are calculated, and the result analysis is obtained.

#### 4.3.2. Model Optimization Process

This model is programmed by LINGO software and solved by GLOBAL algorithm. The model optimization process is shown in Figure 3 below.

First of all, initialization the system, system variables and parameters are defined, input cold, thermal, electric load parameters and *krun*/*kpoll* weight parameters, then, objective function, equipment model constraint and power balance constraint are input.

After that, the LINGO software is used to calculate the optimization model. If there is no feasible solution to the optimization model, then there is no solution to the original optimization problem; otherwise, continue to optimize the process. The results of feasible solutions are compared. If the result of comparison is the minimum of feasible solutions, the global optimal solution is obtained. Instead, continue to compare possible solutions.

Finally, output the optimal result and end the optimization process.

**Figure 3.** Model optimization process.

#### **5. Numerical Simulation and Operation Load Rate Analysis**

#### *5.1. Numerical Simulation and Operation Results Analysis*

This paper selects Jianguo hotel in Zhangjiakou, northern China, to analyze the system's cold, heat and electricity load before summer. At the same time, considering the feasible range of system load rate, the system operating conditions are adjusted by changing the system equipment parameters, thus expanding the system's operational load rate and balancing the system's "source-load" power.

The parameters of electricity, heat and cold load in summer are shown in the following Figures 4–6:

The electricity power load before summer presents obvious "peak and valley" characteristics. The power demand is large during peak hours and low during valley hours. Therefore, it is very suitable for implementing time-of-use electricity price. The specific time-of-use electricity price is shown in Table 1 below.


In summer, the system has a low demand for heat energy, and the load fluctuates. The temperature rises from 10: 00 to 14: 00 noon, and the heat energy load is low. The system has a high demand for hot water and a high heat energy load from 21: 00 to 23: 00 at night.

In summer, the system has a strong demand for cool energy. In the early morning (0:00~5:00), the cool energy load is relatively low, then increases with the temperature rising, reaches the maximum at 14:00, and then continues to decline.

This optimization problem is a mixed integer nonlinear problem. Standard linear method is adopted to solve the multi-objective problem. By changing the weight coefficients of *krun* and *kpoll*, different optimization results are compared, and the global optimal solution is obtained. In terms of optimization model calculation, this model has a total of 1923 variables and 1824 constraints, including 1128 linear variables and 766 linear constraints. In order to ensure that the output result is the global optimal solution, multiple feasible solutions need to be compared, so the calculation time is also relatively long. The calculation time for each group of different weight coefficients needs at least 3–4 h. The optimal results are obtained by changing the *krun* and *kpoll* weight coefficients of multi-objective functions as shown in Figure 7 below:

**Figure 7.** Optimal operation results of different weights in summer.

According to the optimization results of different weights in summer in Figure 7, when the weights (*krun*,*kpoll*) = (0.1,0.9) and (0.2,0.8), the system cannot operate, and the weights of these two points should be discarded. When the weights *krun* = 0.9 and *kpoll* = 0.1, the mixed objective function value *F* is the minimum of 1.0022. At this time, the electric/heat/cool output power of each equipment before summer is shown in the following Figures 8–10:

**Figure 8.** Power output diagram of various equipment before multi-time scale in summer.

As can be seen from Figure 8, the gas internal combustion engine has maintained good electric energy output characteristics. In the 0th to 7th and 23rd hour, the electricity price of the power grid is low, which is suitable for the system to receive a large amount of electric energy of the power grid and meets the requirements of system economy. The gas internal combustion engine is also put into operation, with the starting load rate kept above 55% and the output power stable. In the 8th to 17th hour, the photovoltaic unit is actively connected to the system, the system fully absorbs renewable energy, the power of the power grid is greatly reduced, and the gas internal combustion engine runs stably. At this time, the load rate is kept above 85%, and the electric energy output is stable. During the 18th to 22nd hour of the peak power consumption at night, although the photovoltaic unit will no longer output power, the gas internal combustion engine will remain in a state close to full load output and the power output of the power grid will be limited to a stable range. Electric boilers and electric refrigerators, as units of electric energy consumption, are converted into heat energy and cold energy to supply loads for use. At the same time, the charging and discharging power of the power storage equipment is less throughout the day, which indicates that the system can realize self-absorption, which is beneficial to prolonging the service life of the power storage equipment and the stable operation of the system power.

**Figure 9.** Heat output diagram of various equipment before multi-time scale in summer.

In summer, the demand for heat energy load is relatively low. The heat energy supplied by converting the cylinder liner water directly connected with the gas internal combustion engine into hot water basically bears the base load part of the heat energy load in summer. The electric energy operates twice in one heat energy operation cycle, while the gas internal combustion engine outputs stably in one hour. Therefore, the heat power output by the cylinder liner water is relatively stable in one heat

energy operation cycle. The electric boiler works in the electricity valley period and supplements the heat energy load with the flue gas absorption heat pump. The absorption refrigerator absorbs the surplus heat energy from the cylinder liner water during the summer heat energy valley and converts it into cold energy. The heat storage equipment absorbs and releases less heat energy during one day's operation, which indicates that the system has balanced heat power and stable operation.

**Figure 10.** Cooling output diagram of various equipment before multi-time scale in summer.

In summer, the demand for cooling energy load is large, and the electric refrigerator and absorption refrigerator provide only a small part of the cooling energy load all day long. As the flue gas absorption heat pump directly absorbs the flue gas heat energy of the gas internal combustion engine, it can stably output four times of cooling energy in one cooling energy scheduling cycle. Moreover, the flue gas absorption heat pump is sensitive to load changes, has good load following characteristics, and meets the cooling power load demand of the system.

#### *5.2. Analysis of Operation Load Rate Results*

In terms of system operation load rate, the loads of cold, heat and electricity are independent of each other. Since the system load rate is affected by the loads of electricity, heat and cold at the same time, analyzing the changes of the three to the load rate will directly determine the stable operation of the system. The operation of electric/heat/cold load and load rate is shown in Figures 11–13.

**Figure 11.** Multi-time scale power load-load rate operation.

As shown in Figure 11, the multi-time scale power load and load rate operating points are mostly close to the upper power limit of the decoupled power system, while the load rate operating points close to the upper power limit are close to full load operation, and the system is stable at this time. The closer the operating point is to the upper limit of power, the more thorough the system is running, the greater the utilization rate of power and the higher the economy of the system.

**Figure 12.** Multi-time scale heat load-load rate operation.

The multi-time scale heat load and load rate operation is shown in Figure 12. The summer heat load demand is low and the load fluctuation range is small, so the heat load and load rate operation points are distributed in a scatter line within the thermal energy output range. At the higher load rate, the operating point is close to the lower limit of the heat energy output boundary, but the heat storage equipment has not been put into operation at this time, indicating that the "source-charge" power of the heat energy subsystem is still balanced.

**Figure 13.** Multi-time scale cooling load-load rate operation.

For summer with large demand for cold energy, due to the strong compatibility and output characteristics of the cold energy system, there is still a wide range of load rate operation range selection under the condition of large fluctuation of cold load throughout the day, so that the load rate operation point can be located inside the output range and the cold energy subsystem can operate stably.

To sum up, under the condition that the cold/heat/electric loads are independent of each other and do not interfere with each other before summer, according to the analysis of the load rate operation results of each decoupling subsystem, the "source-load" power of the system is balanced and stable in operation.

#### *5.3. Comparison between Multi-Time Scale and Single-Time Scale Systems*

In order to analyze the difference of time complementarity between systems, this paper compares multi-time scale systems with single time scale systems. In a single-time scale system, according to the principle of minimum adaptability of load scheduling time, i.e., taking the longest scheduling period as the scheduling time, the cold/heat/electric energy systems will be unified into the same time scale scheduling (Δ*t*<sup>1</sup> = 1 h), the multi-time scale electric energy systems will be unified into one hour by 15 min scheduling, the heat energy systems will be unified into one hour by 30 min scheduling, and the original load parameters will be added into one hour to calculate. As the peak power load of a single-time scale system increases nearly 4 times as much as that of a multi-time scale system, and the gas internal combustion engines cannot match the applicable models due to capacity, power and other reasons, it is necessary to increase the number of gas internal combustion engines to 4 in a single-time scale system and keep the load rates of each gas internal combustion engine consistent so as to facilitate systematic analysis and comparison. Then the single time scale system operation results are shown in the following Figures 14–16:

**Figure 15.** Heat output diagram of various equipment before single time scale in summer.

**Figure 16.** Cooling output diagram of various equipment before single time scale in summer.

The power output of each equipment before the single time scale in summer is shown in Figure 14. The gas internal combustion engine operates at a load rate of more than 80% throughout the day, taking up half of the power load throughout the day. The system has completely absorbed the photovoltaic power and greatly reduced the dependence on the grid power when the system is running at noon (the 10th to 15th hour). Electric boilers and electric refrigerators absorb only a small part of electric energy throughout the day. However, power storage equipment is put into operation at certain times to maintain stable operation of the power system.

The heat energy output of each equipment before the single time scale in summer is shown in Figure 15, and the cylinder liner water heat exchanger outputs sufficient thermal energy. Electric boiler and flue gas absorption heat pump output less; The absorption refrigerator only absorbs a small amount of heat energy and converts it into cold energy. During the first 2, 9, 15 and 19, 23 h of system operation, the heat storage equipment absorbs and releases power to maintain the stability of the thermal system.

The cooling energy output of each equipment before the single time scale in summer is shown in Figure 16. The flue gas absorption heat pump provides almost all-day cooling energy load, while the electric refrigerator and absorption refrigerator are only used as auxiliary equipment, thus the cooling energy power of the system is balanced.

In terms of single-time scale system load rate operation, the relationship between system load and operation load rate is shown in the following Figures 17–19:

**Figure 17.** Single time scale power load-load rate operation.

Single-time scale power load-load ratio operation is shown in Figure 17. Although most power load ratio operation points are close to the upper limit of power output of the decoupled power system and the power system operates efficiently, some power load ratio operation points exceed the upper limit of power output boundary. At this time, the power storage equipment is put into operation to absorb some power, otherwise the system cannot operate stably, and at the same time, the energy storage burden of the system is increased.

**Figure 18.** Single time scale heat load-load rate operation.

Single-time scale heat load-load ratio operation is shown in Figure 18. summer heat load is basically linearly distributed. Due to single-time scale adjustment, the heat load of the system is almost twice that of multiple time scales within one hour. At some operating points with high load rate, the lower limit of the power output boundary has been dropped. At this time, the heat storage equipment must be operated to ensure the stability of the system output heat energy.

**Figure 19.** Single time scale cooling load-load rate operation.

The operation of single-time scale cooling load-load ratio is shown in Figure 19. Although the cooling energy scheduling period of single-time scale system is the same as that of multi-time scale system, the system load ratio is determined by the loads of cold, heat and electricity, which makes the overall load ratio of single-time scale system higher and the system is inferior to multi-time scale system in load ratio selectivity.

Comparing the load rates of single-time scale system and multi-time scale system, it is found that the load rates of single-time scale system are kept above 80% throughout the day, while the load rates of multi-time scale system are scattered, ranging from 55% to 100%. On the surface, the single-time scale system load rate is more stable and efficient, but in fact, the single-time scale system maintains a high load rate and the system operates at high load, which will greatly increase the system equipment loss and maintenance cost, which is very unfavorable to the long-term stable operation of the system and also increases the risk of system failure. The load rate of multi-time scale system is stable at 55~70% in low load period, but it can reach full load operation in high load period. The feasible range of system load rate is wider and the system has wider operation selection space.

A comparison of system operation results on multiple time scales and on a single time scale is shown in Table 2 below.


**Table 2.** Comparison table of results for multi-time scale and single-time scale systems.

From Table 2, it can be seen that the results of the multi-time scale system are better than those of the single time scale system in terms of total operating costs and pollutant gas emissions, with a reduction of 15.6% in economic operating costs and 22.3% in pollutant gas emissions. Under the condition of multi-time scale, the system equipment has more flexible time collocation and energy complementary selection, highlighting the multi-energy complementary characteristics of the cold-heat-electricity integrated energy system. However, under the condition of a single time scale, the flexible multi-time scale adjustment mode of the power system and the thermal system is abandoned, and the load burden of the system is increased in the same time dimension, so that the system has to expand the upper limit of equipment capacity or increase the equipment investment when selecting equipment, which is extremely dependent on the energy storage equipment to maintain the "source-load" power balance of the system, thus greatly increasing the investment cost of the system, the difficulty in equipment selection and the maintenance cost of the equipment, and the multi-time scale system has a wider operating load rate range.

#### **6. Summary**

The cold-thermal-electricity integrated energy system takes coupling multiple energy production modes to realize multi-energy complementation, energy step utilization and overall energy utilization rate improvement as the core, thus achieving the purposes of high efficiency, energy conservation, environmental friendliness and compatibility. Because it is close to the energy consumption side, the energy loss in the energy transmission process is greatly reduced. The mutual supplement of various energy sources also reduces the losses between energy production at all levels. At the same time, there are still many difficulties in optimizing the integrated energy system. Based on the multi-objective and multi-time scale optimization problem of the cold-thermal-electricity integrated energy system under the condition of the feasible load rate interval, this paper studies the cold-thermal-electricity integrated energy system topology structure, takes the gas internal combustion engine and the flue gas absorption heat pump as the core, considers the feasible load rate interval of the system under the condition of strong coupling, and adopts the method of decoupling analysis of the system to analyze the cold, heat and electrolysis coupling subsystems. The analysis results show that the selection of the system load rate will directly determine the "source-load" power balance of the system. In order to ensure the stable operation of the system under the condition of independent and uncertain loads, the equipment parameters should be adjusted according to the fluctuation range of each load to ensure the operating conditions of the system, so that the system load rate can fall within the stable load rate range of the cold, heat and electrolytic coupling subsystems. On the issue of multi-objective and multi-time scale optimization, the multi-objective function takes into account the lowest system economic operation cost and the lowest pollutant gas emission. At the same time, considering the

differences of different energy attributes and energy scheduling characteristics, different time scales are selected to model the equipment model and the power model. The power control is sensitive and the scheduling is fast, with 15 min as an operation cycle. The flexibility of heat energy is obvious, taking 30 min as an operation cycle; Cold energy changes slowly, taking 1 hour as an operation cycle. In the aspect of multi-objective problem solving, the unitary linear weighting method is adopted to convert the multi-objective problem into a single-objective problem, and different weighting coefficients are set to optimize the solution, and the above mixed integer nonlinear problem is calculated by LINGO solver. This paper selects a hotel in Zhangjiakou, northern China, for simulation analysis of summer front cold, heat and electricity loads, and compares it with a single time scale system. The results show that the multi-time scale system reduces the economic operation cost by 15.6% and the pollution gas emission by 22.3% compared with the single time scale system. Under the multi-time scale condition, the system equipment has a wider load rate operation range, flexible time allocation and complementary energy selection, and greatly reduces the equipment capacity, investment cost, equipment selection difficulty and equipment maintenance cost of the single time scale system, which is more conducive to the stable operation of the system.

This study uses a static model. When considering multi-time scale optimization, the model parameters change dynamically. By the way, the optimization time of cold energy is long, so it can be considered to reduce the optimized time scale of cold energy appropriately. Although the optimal weight coefficient (*krun*, *kpoll*) is selected after the system optimization results, the results show that when the weight coefficient (*krun*, *kpoll*) is (0.9, 0.1), it performs best, and the target has little impact on the pollution emission control, which is one of the directions that this method can improve. In terms of future work, the system model can be optimized to improve the dynamic parameter changes of equipment, so as to build a real-time coordination and optimization system and improve the stability of comprehensive energy system.

**Author Contributions:** Methodology, B.O., Z.Y. and L.Q.; Writing—Original Draft Preparation B.O.; Guidance and Editing, Z.Y., C.L., L.Q. and D.L.

**Funding:** This work supported by "National Key Research and Development Program of China 2018YFB0905105". **Acknowledgments:** The research gratefully acknowledge the LINGO software to solving optimization. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**





#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **On Minimisation of Earthing System Touch Voltages**

#### **Vaclav Vycital \*, Michal Ptacek, David Topolanek and Petr Toman**

Department of Electrical Power Engineering, Faculty of Electrical Engineering and Communication, Brno University of Technology, 601 90 Brno, Czech Republic; ptacekm@feec.vutbr.cz (M.P.); topolanek@feec.vutbr.cz (D.T.); toman@feec.vutbr.cz (P.T.) **\*** Correspondence: vycital@feec.vutbr.cz; Tel.: +420-541-146-247

Received: 15 September 2019; Accepted: 3 October 2019; Published: 11 October 2019

**Abstract:** Finding cost efficient earthing system design with acceptable level of safety might be quite tedious work. Thus, many earthing system engineers try to find the most suitable design either by employing only their best experience or taking advantage of some more complex optimisation programs. Although both approaches might work well under certain circumstances, they might fail either due to counter-intuitiveness of the specific situation or by misunderstanding of the applied optimisation method, its limitations etc. Thus, in this paper, the earthing system design optimisation problem was addressed by analysing optimisation simulation results together with conducted sensitivity analysis. In the paper, a simple double ring earthing system was optimised while using five different optimisation methods. The earthing system was placed in different horizontally stratified soil models and the earthing system was optimised by minimising touch voltages instead of commonly minimised earth potential rise. The earthing system was modelled by Ansys Maxwell software. Apart from using Ansys Maxwell built-in optimisers, the possible solution space has also been mapped by performing sensitivity analysis with changing the earthing system design dimensions and the results of optimisation were compared and validated. It was found out that the Sequential Non-Linear Programming Optimisation technique was quite superior to the other techniques. Additionally, in most cases, the Ansys Maxwell optimiser was able to found optimal solution; however, in some cases, based on the initial conditions, it might get stuck in local minima or the results might be influenced by the solution noise. Additionally, some quite non intuitive dependencies of earthing system electrodes positions had been found when different spatial dimensions constraints are used.

**Keywords:** earthing; earthing system; optimisation; Ansys Maxwell; sensitivity analysis; touch voltage; non-linear programming

#### **1. Introduction**

The earthing system (ES) is one of the necessary parts of distribution and transmission power networks. The basic requirements on ES are stated in standards EN 50522:2010 [1] or e.g., in its American counterpart IEEE 80 [2]. Apart from other requirements, the designing engineer must ensure that step voltages (SV) and touch voltages (TV) in the vicinity of the ES shall not exceed the safety limits defined by [1], or [2] respectively. Thus, the proposed design by the engineer is always compliant with safety requirements, e.g., permissible touch voltage curve of EN 50522, however it might not be the best possible design for this desired specific site and its constraints. It might be quite common that the initial approach by the engineer is to apply some universal ES design recommendations (e.g., for distribution transformation station use double ring earthing system). The spatial layout of this initial design is usually the same for every such proposed design. Once this initial design is not compliant with the safety requirements, the designing engineer makes some adjustments to this initial design to find a solution that is compliant with the safety limits. However, the final proposed design by the engineer is not optimised within the specific site constraints and there might be still room for

finding better design layout by employing a more complex optimisation approach, i.e., decrease the TV, decrease used material, etc. This designing approach might be especially true in case of the great number of built distribution transformer stations supplying the low voltage distribution networks. In the case of these ESs, the designing engineer very often uses the simplified formulas of an Annex J [1] to get the ES resistance to earth (e.g., ES ring Equation (1))

$$R\_r = \frac{\rho}{\pi^2 D} \ln \frac{2\pi D}{d},\tag{1}$$

where ρ is soil resistivity, *D* is diameter of earthing ring, and *d* is the diameter of ring electrode cross section.

Once the resistance to earth is obtained, the safety assessment is quite straight forward, as the failure current might be obtained from simplified equivalent electrical circuit with the equivalent voltage source method [3]. The ES is then assessed based on the prospective permissible touch voltage that is taken as one half of entire earth potential rise (EPR) [1].

Basically, the abovementioned approach is kind of optimization technique with an objective function

$$\text{minimise } \text{EPR}(\mathbf{x}) \tag{2}$$

Subject to design spatial dimension constraints

$$g\_1 \cdots g\_n \le d\_1 \cdots d\_n \le h\_1 \cdots h\_n \tag{3}$$

where *d*<sup>1</sup> ... *dn* are appropriate ES dimensions and *g*<sup>1</sup> ... *gn*, *h*<sup>1</sup> ... *hn* are the spatial dimension constraints. This is because the designing engineer should mainly focus on decreasing the overall EPR to lower the corresponding touch voltages as much as possible. Decreasing the EPR is most effectively done through decreasing the resistance to earth of the ES design. The described procedure of (2) and (3) in all situations would mean lengthening the ES horizontal dimensions up to the outer boundary of the horizontal constraints and burying the ES near the burial depth lower constraint. This procedure might be also expected, regardless of most often used horizontally stratified soil model if the simplified formula of (1) is used.

Usually, the ES horizontal dimensions are constrained by utility owner limited land. Additionally, burying the ES to greater depth might get costly due to more excavation work needed. Thus, many researchers have been dealing with the ES design optimization problem. Some of the papers were dealing with finding/proposing a method to optimize the design by appropriate ES conductor allocation in the soil with reducing earthing system total EPR only [4]. Other works are aimed on the optimization by minimising the total length of used ES horizontal and vertical conductors [5,6]. Some other researchers also included into the optimization the costs that are connected with the installation process, i.e., the excavation costs as well as the ES conductor material costs [5,7,8]. Another, quite a different approach to ES design optimisation can be found in [9,10]. It is proposed by these papers to optimize the ES by using low resistivity materials. This optimization method is based on the fact that by applying different chemicals near or next to ES rods and strips, the fault current is more effectively driven into the ground due to increased effective surface area of the ES rods. By this method, the ES resistance might be influenced either temporarily or permanently. Lastly, but not used very often might be the optimisation by finding optimal distance to some other external grounding grid [11].

Generally, two things are necessary for the optimization procedure. The first thing is some objective function that is used to quantify if the so-called optimum has been reached. The second necessary thing is some kind of step direction function that determines the direction towards the so-called optimum. The usage of some different objective functions has been described in the previous paragraph. However, in the case of the step direction function, the aforementioned papers [4–11] mainly use two different approaches or their combination. The first approach is by performing parametric

analysis (i.e., sensitivity analysis). This approach is based on the fact that, if there are only expected smooth changes in the solution, by choosing appropriate parametric step, the whole solution space can be mapped and this solution space is later searched for the so-called optimal solution according to the objective function. The results of this parametric analysis can be also used as a recommendation for the designing engineer if properly depicted [12]. The second approach that was identified in the reference is the use of some more complex optimization approach e.g., using non-linear programming optimization methods, like genetic algorithm optimization, particle swarm optimization, and hybrid particle swarm genetic algorithm optimization. In general, ES design optimization is a spatial problem where the use of the mentioned methods can decrease the overall time that is needed for the parametric analysis. In some cases, performing parametric analysis with fine steps might get quite tedious and time consuming.

Quite a novel approach to earthing system optimization problem might be the optimization through Quantified Risk Analysis [13–15]. The objective function of this approach is the quantified real risk (i.e., the resulting risk of death consisting of the probability of hazardous TVs level and the coincidence probability [13]). This approach is mostly assumed as an alternative to the conventional 'worst case' earthing system approach [1], which might be more strict or more loose, depending on the situation. However, the authors are yet unaware of publication using this approach together with some kind of step direction function and, also, the so-called 'optimum' for this method might be questionable due to diverse risk perception and acceptable risk level laws.

Not many papers have been dealing with the earthing system optimisation through the minimisation of TVs in case of horizontally stratified soil model. For example, in paper [4] the total EPR together with the TV was analysed. The 'compression ratio' function was proposed as a step direction function. The proposed approach is best suited for horizontal earthing mats with parallel and perpendicular conductors and it can be easily used for conducting parametric analysis. In the case of searching for optimal compression ratio, a more sophisticated step direction function different from simple first or second order gradient based approach might be necessary. Additionally, only the horizontal separation of parallel earthing conductors is optimised without any mean of incorporating earthing mat or separate electrodes burying depth.

Thus, in this paper, an optimisation analysis is introduced for the allocation of separate earthing system electrodes in earth through using Ansys Maxwell Optimetrics Component [16]. The Optimetrics Component allows for the user to search for ES design optimal dimensions. The prospective TV was assumed as the source of real associated risk with the earthing system design and it was chosen as the objective function of the optimisation. A comprehensive parametric/sensitivity analysis was also carried out with changing some of the model parameters to assess the correctness of obtained results by the built-in optimisers in Ansys Maxwell [16]. By this approach, the solution space was searched first and the results that were obtained from sensitivity analysis were compared with results from the built-in Maxwell optimisers.

#### **2. Ansys Maxwell Optimisers**

Ansys Maxwell [16] is a software module incorporated in the whole simulation software ANSYS. Ansys Maxwell is software for the simulation of low frequency quasi stationary electromagnetic fields (EMF). The EMF can be determined for general spatial electromagnetic problems as well as for only two-dimensional problems. The solution of Maxwell equations is managed through the finite element method. Basically, the user has to define/import the required design model, define sources, boundary and meshing condition, and the solution is determined by approximating the EMF distribution through the mesh of finite elements [17]. Additionally, with the use of the built-in field calculator a user defined quantities might be obtained from postprocessing of solved EMF data.

Ansys Maxwell presents an Optimetrics component that enables the user to perform parametric/sensitivity/sweep analysis. Furthermore, a component for optimisation is available. For the Optimetrics analysis it is necessary to define *user defined variables*, i.e., model length and

any other dimensions, and *user defined EMF quantities*, like electric potentials at certain points, etc. The user defined variables are later used in the optimisation analysis as the parameters that are changed in iterative manner to find optimum design. The user defined EMF quantities are used for definition of objective function of the optimisation. The software features six different options for model optimisation [16]:

#### *2.1. Quasi Newton (QN)*

The Quasi Newton optimizer is a gradient based optimizer numerically determining the derivatives of quadratically approximated higher order functions of actual local objective function dependency. This optimizer suffers from the possibility of stucking in local minima, and due to finite element method nature of the results also stucking due to high level of solution noise. As the Hessian matrix increase rapidly with the number of optimised design variables the computational time might increase rapidly and thus this optimizer is recommended for optimisation of 1–2 design variables only.

#### *2.2. Pattern Search (PS)*

The Pattern Search optimizer uses the grid based simplex search with simplices in the form of tetrahedrons. It is non-gradient based optimisation search, thus it is less dependent on the solution noise. However, it might take longer to find the optimum. The use of only three optimised design variables is recommended.

#### *2.3. Sequential Non-Linear Programming (SNLP)*

The SNLP optimizer is more complex one. Basically, the optimizer creates a response surface by using Taylor series approximation of finite element analysis. By this approximation, the optimizer can find locations of improving points and determine next step direction. The SNLP optimizer is supposed to be more accurate and reliable because it is not constrained by the problem of stucking in local minima like QN. That is also caused by the possibility to overcome the local minima by taking larger steps within the optimisation variable constraints limits.

#### *2.4. Sequential Mixed Integer Non-Linear Programming (SMINLP)*

SMINLP optimizer is basically same as SNLP. However, this optimizer allows for mixing discrete and continuous optimised design variables.

#### *2.5. Genetic Algorithm (GA)*

The GA optimizer is stochastic optimizer. It runs many iterations with random selection of next searched points. It uses an iterative process with initial generation/parents/children and mutated population. The better performing individuals are chosen for the generation of the next generation in each iteration. Through this random search process, the solution space is searched in a structured manner and the optimum is found. However, the disadvantage is also that the non-improving designs are searched in the random selection process and thus this approach takes many more iterations and is thus slow.

#### *2.6. Matlab Optimiser (MO)*

Ansys Maxwell also facilitates the possibility to pass parameters to Matlab software (Matlab 2016b, MathWorks, Natick, USA) and invokes Matlab built-in optimisation functions and uses them to optimize the Ansys Maxwell design.

#### **3. Earthing System Optimised Model**

The feasibility of using Ansys Maxwell software for earthing system design optimisation was assessed on simple double ring earthing system, as in Figure 1. The first inner ring with a diameter *D*<sup>1</sup> was buried at depth *h*<sup>1</sup> and the outer ring with diameter *D*<sup>2</sup> was buried at depth *h*2. The earthing system was made of FeZn strip with cross section 30 × 4 mm2. The earthing system was modelled in Ansoft Maxwell [16,17] with horizontally stratified soil model with two layers—the surface one with resistivity ρ<sup>1</sup> and thickness *H* and the bottom layer with resistivity ρ2. The soil was modelled by semi-sphere with radius of 200 m and the analysis percent error was set to 1%. For the illustrative purpose, the earthing system model in Figure 1a,c is depicted with an improper soil scale. The earthing system was excited by a DC current [17] of 30 A (red arrow in bottom right picture of Figure 1c) injected at the centre of the ES, where the ES conductor is brought up above the ground and it is considered as accessible for touch.

**Figure 1.** Earthing system model, spatial view (**a**), top view (**b**), and side view (**c**).

#### **4. Earthing System Sensitivity Analysis Results**

A comprehensive sensitivity analysis was carried out before the optimisation of ES design by Ansys Maxwell. In this analysis, the ES was placed in seven different soil models, including 'High on Low' (HoL), Uniform, and 'Low on High' (LoH) soil models, as in Table 1.


**Table 1.** Sensitivity analysis soil models.

For each of used soil model in the sensitivity analysis, the diameters of the ES were changed as inner ring *D*<sup>1</sup> = 1, 2, and 3 m, and outer ring always greater than inner ring *D*<sup>2</sup> = 4 and 5 m. In the case of burial depths of both rings *h*<sup>1</sup> and *h*2, seven burial depths were used for each ring as 0.3–1.5 m with step of 0.2 m. I.e., for one set combination of ES rings diameters *D*<sup>1</sup> and *D*<sup>2</sup> all 49 combinations of ES design burial depths were modelled by Ansys Maxwell and the results were used in the sensitivity analysis. By this way, about 2000 ES design variations were analysed in this analysis. For each of these ES design variations, the potential profile on the earth surface along x direction (Figure 1b) was read and the TV distribution have been determined in the vicinity of the ES. Although the potential profile was read with step of 0.25 m, for the purpose of optimisation evaluation throughout this paper four main performance values have been determined—*TV*1m, *TV*2m and *TV*3m as TV 1, 2, and 3 m apart

from the centre of the earthing system, respectively (i.e., potential difference between EPR and point 1, 2 and 3 m apart). Additionally, the total EPR was read. For set diameter *D*<sup>2</sup> (used here to represent the dimensional constraint), the optimum burial depth was manually searched whilst changing other three dimensions *D*1, *h*1, and *h*2. Incorporating the optimisation of inner ring diameter *D*<sup>1</sup> whilst manually searching for optimum design is quite challenging as the number of analysed variations increase rapidly. The results of the sensitivity analysis are introduced in Figures 2–5 and Tables 2–5.

**Figure 2.** TV1m solution surfaces for ES with diameters *D*<sup>1</sup> = 1–3 m, *D*<sup>2</sup> = 4 m, ρ<sup>1</sup> = 500 Ωm, ρ<sup>2</sup> = 100 Ωm, *H* = 2 m.

**Figure 3.** Touch voltage (TV) surfaces, 500 Ωm/100 Ωm/1 m, *D*<sup>2</sup> =4 m, (**a**) *TV*1m surface; (**b**) *TV*2m surface.

**Figure 4.** EPR, 500 Ωm/ 100 Ωm/ 2 m, *D*<sup>1</sup> = 2 m, *D*<sup>2</sup> = 4 and 5 m.

**Figure 5.** TV profiles along x direction for HoL model ρ<sup>1</sup> = 500 Ωm, ρ<sup>2</sup> = 100 Ωm, *D*<sup>1</sup> = 2 m, *D*<sup>2</sup> = 4 m, *h*<sup>1</sup> = 0.3 m; (**a**) surface layer thickness *H* = 2 m; and, (**b**) Surface layer thickness *H* = 1 m.


**Table 2.** Sensitivity analysis results for ρ<sup>1</sup> = 500 Ωm, ρ<sup>2</sup> = 100 Ωm, *H* = 2 m, *D*<sup>2</sup> = 4 m, *D*<sup>1</sup> = 2 m.

**Table 3.** Sensitivity analysis results for ρ<sup>1</sup> = 500 Ωm, ρ<sup>2</sup> = 100 Ωm, *H* = 2 m, *D*<sup>2</sup> = 4 m, *D*<sup>1</sup> = 3 m.



**Table 4.** Sensitivity analysis results for ρ<sup>1</sup> = 100 Ωm, ρ<sup>2</sup> = 500 Ωm, *H* = 2 m, *D*<sup>2</sup> = 4 m, *D*<sup>1</sup> = 2 m.

\* Design with *h* = 1.5 m EPR selected.

**Table 5.** Sensitivity analysis results for ρ<sup>1</sup> = 500 Ωm, ρ<sup>2</sup> = 100 Ωm, *H* = 2 m, *D*<sup>2</sup> = 5 m, *D*<sup>1</sup> = 2 m.


In Figure 2, the general overview of the results of sensitivity analysis is depicted. For set dimensional constraint *D*<sup>2</sup> = 4 m, three *TV*1m surfaces are depicted for different inner ring diameter *D*<sup>1</sup> and each surface is consisting of results of all 49 burial depth variations. In this compact form, the optimum burial depth of both rings can be immediately read. Within the dimensional constraints of *h*1, *h*<sup>2</sup> (being in range of 0.3–1.5 m) and *D*<sup>2</sup> (being fixed as 4 m), the 'global' minimum of the solution surface can be expected for *D*<sup>1</sup> = 2 m, *h*<sup>1</sup> = 0.3 m, and *h*<sup>2</sup> = 1.5 m. In case of *D*<sup>1</sup> = 3 m, a remarkable saddle point (*h*<sup>1</sup> = 0.7 m, *h*<sup>2</sup> = 0.6 m) is formed by a significant change in slope direction, where two areas with local (*h*<sup>1</sup> = 1.5, *h*<sup>2</sup> = 0.5 m) and global (*h*<sup>1</sup> = 0.3, *h*<sup>2</sup> = 1.5 m) minima are formed. Although this saddle point might not be as much remarkable for all simulated inner ring diameters, it is still present among all the surfaces.

Similar TV surfaces were also obtained for other objective functions in form of *TV*2m and *TV*3m. The results of these analysis are summarized in tabular form in Tables 2–5. In this compact tabular form, only the remarkable points of the sensitivity analysis surfaces are listed for different soil models and for different ES dimensions whilst preserving the clarity of the results. The burial depths of ES designs were selected for two extreme cases and for all three analysed *TV*xm objective functions:


In this way, six different ES designs for set dimensional constraints of *D*1, *D*2, and selected soil model are listed and they can be compared for different objective functions *TV*xm, EPR objectives etc.

In Tables 2 and 3, a slight difference in TVs for HoL soil model might be observed when changing the inner ring diameter *D*<sup>1</sup> from 2 m to 3 m and with set constraint of outer ring *D*<sup>2</sup> = 4 m. Even though the best possible design assessed based on EPR would be with both rings in the bottom layer *h*<sup>1</sup> = *h*<sup>2</sup> = 1.5 m and *D*<sup>1</sup> = 2 m, based on *TV*1m the best design would be a different one with only the outer ring in the bottom layer and the inner ring in the surface layer *h*<sup>1</sup> = 0.3 m, *h*<sup>2</sup> = 1.5 m and *D*<sup>1</sup> = 2 m. Additionally, if *TV*2m and *TV*3m objectives would be assessed a completely different designs are optimal. These findings might not be as much intuitive. However, with the increasing number of possible assessed option, the situation might get even more tricky.

For example, if additional constraints (i.e., different objective function or less favourable soil model) would be solved, a kind of unexpectable design solution might have been found. This can be apparent from Figure 3a,b, where the TV surface of HoL soil model with surface layer thickness of only 1 m is depicted for objectives *TV*1m and *TV*2m and for constraint of outer ring diameter *D*<sup>2</sup> = 4 m. As it

might be expected by some experienced engineers, the better designs are with an outer ring burial depth greater than the thickness of the surface layer. The overall best design, as already stated in the previous paragraph, is with *D*<sup>1</sup> = 2 m, *h*<sup>1</sup> = 0.3 m, and *h*<sup>2</sup> = 1.5 m. However, this might not be true if other constraints are applied, e.g., dimensional constraint of both rings maximum burial depth of only 0.9 m and both *TV*1m and *TV*2m should be assessed. In the region of constraint burial depth of only 0.9 m, there is an evident change of pattern in optimal burial depths (between (a) and (b) of Figure 3) where the optimum design depending on the objective function is either for inner ring close to the surface with *h*<sup>1</sup> = 0.3 m and *h*<sup>2</sup> close to bottom layer, or the inner ring close to bottom layer and outer ring close to the surface, respectively. Additionally, design with increased inner ring diameter of 3 m gives better results than 2 m as in the case of *TV*1m objective. From both figures it is evident that if both *TV*1m and *TV*2m objectives need to be evaluated, the contradiction of the solution will need some sort of more complicated addressment, e.g., weighting of risk, etc.

Additionally, based on the results of depicted sensitivity surfaces, as in Figure 3a, the positive effect of earthing electrode grading can be assessed. For example, for the analysed arrangement and objective of *TV*1m, burying the inner ring only to a shallow depth of 0.3 m led to decreasing the risk. The rest of the designs are either comparable, or even worse. Burying both rings to the bottom layer is worse by slightly more than 40% than the optimal shallow inner ring placement (*h*<sup>1</sup> = 0.3 and *h*<sup>2</sup> = 1.5 m). The TV surface for inner ring diameter of 1 m for HoL has yielded in almost all cases to higher TVs and is thus excluded here.

The Uniform and LoH soil models are almost excluded from presented results, as the optimisation of ES design through changing design dimensions had greater observable dependency in TVs and EPRs for HoL soil model. Table 4 presents the results for LoH soil model. In case of LoH and Uniform soil models the TV differences between the min and max TV options are of small values and thus no big difference was found while searching for optimal solution. Additionally, the EPRs of different designs are almost identical. The abovementioned findings are true if the outer ring diameter is kept constant as a dimensional constraint. If that diameter would be even larger, the differences in TVs would smoothen even more. It is worth pointing out that, in the case of Uniform soil model, the results are even more convenient than in the case of the LoH model.

Increasing the outer ring diameter would of course lead to more suitable ES design, as might be apparent from Table 5, in comparison to Tables 2 and 3, where design with increased outer ring diameter of 5 m is presented (compared to 4 m in Tables 2 and 3). To complete the ES behaviour overview, Figure 4 depicts the EPR surface for different burial depths and outer ring diameters.

Figure 5a,b depict a kind of complementary and expanding overview of ES behaviour. In these figures the potential profiles along the x direction (Figure 1b) on the surface is depicted for ES designs with changing only the outer ring burial depth while keeping the rest of the dimensions constant. From the figures the *TV*1m, *TV*2m, and *TV*3m can be obtained. From Figure 5a, it is obvious that the optimum design in case of objective *TV*1m is with deep outer ring *h*<sup>2</sup> = 1.5 m and with shallow inner ring *h*<sup>1</sup> = 0.3 m. However, this chosen design is by no means also optimal when also assessed to *TV*2m.

From Figure 5b the positive effect of burying the outer ring into the bottom layer in HoL soil model with *H* = 1 m is observable, where the TV profile had dramatically flattened.

#### **5. Ansys Maxwell Earthing System Optimisation Results**

The optimisation of the ES design was also carried out by Ansys Maxwell Optimetrics component for soil models number 1, 3, 4, and 5 (Table 1). The ES was always modelled with some 'initial' dimensions, denoted by the subscript '*i*' (i.e., the *D*1, *h*1, and *h*<sup>2</sup> as *D*1i, *h*1i and *h*2i) and through the optimisation the 'optimised' dimensions have been obtained and are denoted by the subscript '*o*' (i.e., *D*1o, *h*1o, *h*2o). The outer ring diameter *D*<sup>2</sup> was set as a dimensional constraint equal to 4 m in most of the simulations. The burial depth constraint was set as in region 0.3–1.5 m for both rings and the maximum inner ring diameter constraint was set as *D*1max = 0.9x*D*2, thus so the inner ring is always smaller than the outer ring. The minimization of *TV*1m was set as the optimisation objective. The optimisation was

carried out by five different optimisers—QN, SNLP, SMINLP, PS, GA. After the optimisation process some other parameters have been also recorded for possible comparison and benchmarking of different optimisers. I.e., *TV*2m, *TV*3m, the total *EPR* of optimised ES design, total time needed for finding the optimal design *t*, the number of performed simulations in total *ItN* (optimisation iterations), and a number of iterations (out from total of *ItN* iterations) where the optimal design was found *ItO*. As an optimisation analysis stopping criterion was also set the maximum number of optimisation iteration that was in most cases 300, or in some cases only 100. The optimisation results, together with the performed stopping option, are listed in Table 6 (Parts 1 and 2). The stopping option is in the last column 'Status' as either one of the following three options:



*Energies* **2019** , *12*, 3838

**Table 6. Part 1.** Ansys Maxwell

optimisation

 analysis results. **Part 2.** Ansys Maxwell

**Part 1.**

optimisation

 analysis results.



100 *onlylight orange . \* Uniform and LoH soil modelslight blue*

.

**Part 2.**

For HoL soil model 1, it was expected that the optimum burial depth for both rings with fixed outer ring diameter *D*<sup>2</sup> = 4 m (Figure 2) should be 0.3 and 1.5 m for *h*<sup>1</sup> and *h*2, respectively, and inner ring *D*<sup>1</sup> = 2 m (as was found in sensitivity analysis Table 2, Table 3 minimum *TV*1m). For simplicity and comparability, first two optimisation analysis had been performed just with only burial depths *h*<sup>1</sup> and *h*<sup>2</sup> being optimised (lines 26, 27 of Table 6). This is also due to the fact that the simulated designs in sensitivity analysis were simulated with quite great step of inner ring diameter of 1 m, as the possible number of designs rapidly increase when a lower step would had been used. These two optimisation analysis were performed with different initial conditions. In both cases, Ansys Maxwell found ES design close to the expected optimal solution. In the second case, the outer ring depth was only 1.4 m where the optimiser evaluated visiting the depth of 1.5 m as unnecessary. This might be caused by the level of solution noise, i.e., with the set analysis percent error of 1% (Section 3) the changes in solution data of 5–15 V might be expectable due to randomness of creation of finite element mesh and also due to not sufficiently fine mesh elements. Better results might be expected for percent error of about 0.3% [17], however this will cause about three times or even greater increase in solution time.

In the rest of the simulation, the inner ring diameter was also optimised. As the sensitivity analysis was only done with 1 m step of the inner ring diameter, the results of optimisation cannot be really matched to sensitivity analysis results. However, from the comparison of *TV*1m surfaces from Figure 2 it can be expected that the optimal design should has inner ring diameter in 2–3 m region, and the burial depths as *h*<sup>1</sup> = 0.3 m and *h*<sup>2</sup> = 1.5 m. This expectation had been met in most of the results for SNLP optimiser. However, in case of the other optimisers, kind of invalid and different results had been found. This might have happened in some cases due to stucking in actual local minima, stucking in local minima caused by the solution noise, failing due to Ansys Maxwell Adaptive Mesher, etc.

In the case of GA optimiser, more extensive analysis of correct setting of the optimiser might lead to better results. For example, out of four different settings of the optimiser, the optimal design was found in only one (line 10). However, this optimisation had been stopped by user after more than two days of continuous simulation. From the optimisation results it had been found that the optimiser visited the location of global minimum (*D*<sup>1</sup> ~ 2 m, *h*<sup>1</sup> ~ 0.3 m, *h*<sup>2</sup> ~ 1.5 m) only about three times out of total 1134 optimisation iterations. The optimal design was already found after about the first half of all iterations in the 664th iteration, however the optimiser did not evaluate it as the global optimum and was further searching for better solution.

The QN and SMINLP optimisers only found a local minima or did not find the expected optimum design at all in their all 8 performed simulation. Even though the suggested results by the optimisers are quite close to the global optimum, the global optimum still was not reached. The PS optimiser reached the global optimum in only one case out of three and the results seemed to be heavily reliable on the design initial condition. In case of all these three QN, PS, and SMINLP optimisers, there are not many options in Ansys Maxwell software that the user can change to improve the performance of the optimiser. One possible way might be to decrease the analysis percent error, which might help to overcome the optimiser problem of stucking in local minima also caused by solution noise. However, a corresponding increase in solution time is expectable.

The SNLP optimiser tended to give more reliable results, so a more extensive analysis was conducted for this optimiser only focusing mainly on the optimiser immunity to setting different initial conditions of the simulation. An outer ring diameter constraint of 4 m was set in most of the simulations. It was found out that the local minimum becomes significant once the inner ring diameter is close to 3.5 m when a saddle point (Figure 2, comments in Section 4) separates it from the global minimum. This assumption was also clarified by performing additional sensitivity analysis with inner ring diameters of 1.5, 2.5 and 3.5 m (Figure 6 solution surface for *D*<sup>1</sup> = 3.5 m). However, in case of this local minimum the *TV*1m value is about 15–25 V higher that the expected global minimum in the region of *D*<sup>1</sup> = 2–2.5 m. Thus, from results of lines 24 and 28 in (and also 2 QN and 5 PS) Table 6, it can be seen that, in the case of unfavourable initial conditions, the global minimum is not always found and the optimiser might be unable to overcome this problem on its own.

**Figure 6.** TV1m solution surfaces for ES: 500 Ωm/100 Ωm/2m, *D*<sup>1</sup> = 3.5 m, *D*<sup>2</sup> = 4 m.

Lastly, in most cases, the maximum number of 300 optimisation iterations was enough to find the optimal design. In the case of SNLP optimiser, usually slightly more than 100 iterations were necessary to find the global minimum of the optimisation problem. However, if for example the maximum number was set only to 100 iterations, only a solution close to global optimum was found and the 100 iterations limit might be too restricting, i.e., in case of line 16.

#### **6. Conclusions**

In this paper, the appropriateness of using the Ansys Maxwell optimizer to optimize ES design was demonstrated on couple of examples. The use of the SNLP optimiser might be assessed as satisfactory and superior to other Ansys Maxwell optimisers. In most cases, the SNLP method found the optimum design, however in the case of unfavourable initial conditions it is still prone to stucking in local minima. This problem can be overcome by conducting the optimisation analysis for different initial conditions to see whether the found result is kind of reliable. In the case of other compared optimisers, different settings might lead to more reliable results, however the increase in solution time as well as more experienced engineer (user) might be necessary.

Another novel approach that is presented in this paper is the optimisation of the design by minimising touch voltages instead of most often used EPR. Unlike simple EPR minimisation would lead to burying of ES in greater depth it was pointed out that this approach might not necessarily lead to optimal layout also assessed based on actual touch voltages. As the real risk that is associated with the ES is driven by the touch voltages instead of EPR, using of touch voltages is more precise than simple optimisation of EPR. However, when the design is optimised to minimise TVs, it was found out that the optimum ES design might get quite counter intuitive considering the wide variety of options, like different dimensional constraints, different soil models, or even different objective functions ~ transferred voltages, step voltages. Thus, as was presented in this paper, the use of software with more complex optimisation methods (non-linear programming, like SNLP method) might be beneficial, if not even necessary.

**Author Contributions:** Conceptualization, V.V.; Funding acquisition, P.T.; Investigation, V.V.; Supervision, D.T. and P.T.; Validation, M.P. and D.T.; Writing—original draft, V.V.; Writing—review & editing, M.P., D.T. and P.T.

**Funding:** Authors gratefully acknowledge financial support from the Technology Agency of the Czech Republic (project No. TN01000007).

**Acknowledgments:** Authors gratefully acknowledge the Centre for Research and Utilization of Renewable Energy (CVVOZE) where this research work was carried out.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Using the Thermal Inertia of Transmission Lines for Coping with Post-Contingency Overflows**

#### **Xiansi Lou, Wei Chen and Chuangxin Guo \***

College of Electrical Engineering, Zhejiang University, Hangzhou 310027, China; louxiansi@zju.edu.cn (X.L.); chenwei\_ee@zju.edu.cn (W.C.)

**\*** Correspondence: guochuangxin@zju.edu.cn

Received: 31 October 2019; Accepted: 17 December 2019; Published: 20 December 2019

**Abstract:** For the corrective security-constrained optimal power flow (OPF) model, there exists a post-contingency stage due to the time delay of corrective measures. Line overflows in this stage may cause cascading failures. This paper proposes that the thermal inertia of transmission lines can be used to cope with post-contingency overflows. An enhanced security-constrained OPF model is established and line dynamic thermal behaviors are quantified. The post-contingency stage is divided into a response substage and a ramping substage and the highest temperatures are limited by thermal rating constraints. A solving strategy based on Benders decomposition is proposed to solve the established model. The original problem is decomposed into a master problem for preventive control and two subproblems for corrective control feasibility check and line thermal rating check. In each iteration, Benders cuts are generated for infeasible contingencies and returned into the master problem for adjusting the generation plan. Because the highest temperature function is implicit, an equivalent time method is presented to calculate its partial derivative in Benders cuts. The proposed model and approaches are validated on three test systems. Results show that the operation security is improved with a slight increase in total generation cost.

**Keywords:** security-constrained OPF; preventive control; corrective control; post-contingency overflow; transmission line; thermal inertia

#### **1. Introduction**

The security-constrained optimal power flow (SCOPF) is an essential tool for making a day-ahead and real-time generation plan [1]. From the perspective of control measures, the SCOPF model can be divided into the preventive SCOPF (PSCOPF) model and the corrective SCOPF (CSCOPF) model. The former model is widely used in current operations of large-scale [2] and island power systems [3]. The operating point obtained by this model can guarantee the safe operation of systems under all credible contingencies. How to filter the contingency set determines to a great extent the performance of PSCOPF model [4]. For given contingencies, results of the PSCOPF model are conservative with the high generation cost. The risk conception has been introduced in [5] to enhance its economics. In [6], an identification method of superfluous constraints has been proposed and the model of PSCOPF can be simplified and efficiently solved. With the large-scale access of sustainable energy and frequent occurrences of natural disasters, the operation environment of power systems becomes more complex and the feasible region of PSCOPF model may not exist in some difficult operating conditions.

For seeking the lower operating cost and more flexible dispatch plans, the CSCOPF model has been put forward. In this model, corrective measures such as the unit rescheduling and load shedding are utilized after the contingency occurrence to maintain the transmission security of systems [7]. There have been some trails about using the CSCOPF model in energy management systems to help operators make timely and reasonable dispatch decisions [8]. However, two main problems hinder its wide application. The first problem is the heavy computational burden caused by the

huge number of contingencies especially for large-scale systems and coupling constraints between the preventive control and corrective control. Many efforts have been devoted to shortening the solving time. In [9], Benders decomposition has been applied to solve the CSCOPF model, and the computational complexity has been analyzed. Results showed that the computation speed was significantly improved without sacrificing the accuracy, whether in the serial or parallel computing environment. In [10], an iterative approach that comprises four modules has been proposed, and its performance was better than the direct approach and Benders decomposition. A hybrid method has been used to solve the CSCOPF model in [11] where the maximum feasible region was randomly searched through evolutionary algorithms and then deterministic solutions were provided by the interior-point method. A similar solving strategy has also been adopted in [12], and the optimal coordination of the preventive control and the corrective control was obtained.

A more significant problem of the CSCOPF model is the insecurity of power systems during the post-contingency stage. After the contingency such as line failures or generator outages, conventional corrective measures like the unit rescheduling cannot be immediately implemented due to the limit of ramping rates. The system still operates under the preventive control while violations of line power flow and bus voltages may happen. Once cascading failures are triggered, the effectiveness of pre-made corrective plan is influenced, and the scope of the original accident could be extended. Some previous work has been done to improve the security of the CSCOPF model. An optimal locating method for support generator units has been proposed in [13] to improve ramping abilities and enhance the robustness of systems. In [14], quick-start units have been utilized in corrective actions for reducing the duration of post-contingency stage. References [15,16] have discussed that the fast-response distributed battery energy storage can be used to alleviate post-contingency overloads and reduce the power flow below the short-term emergency rating. In [17], the state of charge has been considered in the distributed energy storage model and results showed that the generation cost increased as the initial charging state declined. The post-contingency demand response has been used to relieve overflows incurred by the renewable generation fluctuation and "N−1"contingencies in [18]. References [19,20] have analyzed the risk of cascading failures caused by post-contingency overloads in the alternating current (AC) and direct current (DC) power transmission network, while the thyristor controlled series capacitor and the multi-terminal direct current have been applied to minimize the load shedding and generation rescheduling, respectively. In addition, the model predictive control (MPC) provides a higher perspective beyond static and open-loop approaches. It is a class of strategies that utilize a process model to establish a control sequence for controlling the future behavior of systems over a horizon [21]. Under this framework, the generation redispatch, load shedding, and regulating transformers have been applied to alleviate emergency thermal loads in [22,23]. Moreover, the energy storage and curtailment of renewables have also been identified as corrective actions in [24].

In this paper, the DC power flow is utilized while post-contingency line overloads are the insecurity problem discussed and solved. The essence of this problem can be described as the lack of short-term transfer capacities of power grids. Compared with approaches such as using quick-start generators, distributed battery energy storage, demand response, and other power flow controllers, the method of improving line capacities could be more direct and effective. The line thermal rating provides a novel perspective instead of conventional power flow limits. The static power flow rating obtained under the worst environment condition is usually conservative. Reference [25] has provided a report of two pioneering schemes in the U.S. and the U.K. where the real-time thermal rating was applied. A new power flow formulation considered the line electro-thermal coupling effect has been introduced in [26], and the error in power transfer evaluation decreased by about 20% by using this method. In [27–29], the thermal rating of transmission lines has been utilized to increase the threshold value of wind power integration and avoid the unnecessary tripping of pre-selected generation assets. Besides the static rating difference, transmission lines have the thermal inertia and corresponding time constants range between several mins and ten mins [30]. This means that dynamic variations of conductor temperatures cannot be ignored especially for the short-term post-contingency stage

discussed in this paper. Reference [31] has demonstrated that the neglect of transient thermal behaviors of lines caused the underestimation of power transfer capabilities and led to misoperations. In [32,33], the concept of electrothermal coordination has been proposed and the numerical analysis through a number of case studies validated its benefits in augmenting power transfer capability, emergency control, congestion management, and improving the system loadability. A temperature dependent power flow model has been established in [34] where the dependence of line resistances on conductor temperature was taken into account. In [35], the dynamic electrothermal effect has been brought into safety constraints of the post-contingency power flow control to excavate the overload endurance capability of lines.

When the dynamic thermal process of lines is considered in optimal power flow models, the heat balance equation (HBE) which describes variations of conductor temperatures must be integrated. Because the radiated heat loss rate is proportional to the fourth power of the conductor temperature, the HBE is a nonlinear differential equation. How to solve those electrothermal coupling optimization problems becomes a challenge. In [32,33], a modified Euler method has been used to discretize the HBE and then the problem was broken down into several linearized subproblems whose solutions were iteratively refined and linked. An approximate quadratic relationship between the temperature and the square of current has been derived by some reasonable approximations in [34], and the analytical solution of HBE was obtained for a step change in currents. Reference [35] has proposed that, through transforming the HBE into algebraic difference equations by a numerical integration method, the optimal solution can be realized by combining those transformed equations with other algebraic equations.

This paper proposes that the thermal inertia of transmission lines can be used to cope with post-contingency overflows in the conventional CSCOPF model. According to the system behavior, the post-contingency stage is divided into a response substage and a ramping substage in time sequence. In the previous substage, corrective measures have not been implemented while the power flow shows a step change. In the latter substage, the power flow variation is approximately described by a linear function with the adjustment of unit active power outputs. Based on the assumption of constant resistances and the linearization of the radiated heat loss function, analytical solutions about real-time temperatures in both substages are obtained from the HBE. An enhanced security-constrained OPF (ESCOPF) model is established on the basis of CSCOPF model where thermal rating constraints are integrated in the post-contingency stage. The objective of the established model is to minimize the generation cost. In order to effectively deal with this multi-stage optimization problem containing the quadratic objective function and nonlinear constraints, a solving strategy based on Benders decomposition is proposed. The original problem is decomposed into a master problem of the preventive control and two subproblems of the corrective feasibility check and the thermal rating check. Due to the independence between each contingency, the parallel algorithm can be implemented in solving those two subproblems. An equivalent time method is presented to build an explicit relationship between the highest conductor temperature and unit power outputs. In each iteration, corresponding Benders cuts are generated for infeasible contingencies and returned into the master problem. The final generation plan is obtained until those two checks are satisfied for all contingency scenarios. Numerical simulation results validate that the system security is markedly improved by considering the line thermal inertia in the post-contingency stage and the generation cost is just slightly higher than the result of CSCOPF model.

The main contributions of this paper are summarized as follows:

(a) Two typical temperature variations of post-contingency overload lines are analyzed and their analytical solutions are obtained. If the post-contingency power flow is significantly higher than the value in the corrective stage, the stationary point of the conductor temperature exists. Otherwise, its variation is monotonous.

(b) An ESCOPF model is established to minimize the system generation cost. Thermal rating constraints for line conductors are integrated in the post-contingency stage to avoid cascading failures. In addition, then, the unit rescheduling is implemented as the corrective measure to remove long-term overloads on transmission lines.

(c) A Benders decomposition based solving strategy is proposed to solve the ESCOPF model, which is divided into a master problem and two subproblems. In order to derive the partial derivative of the highest conductor temperature with respect to the unit active power output, an equivalent time method is presented.

The remainder of this paper is organized as follows. Section 2 analyzes temperature variations of post-contingency overload lines and solutions of HBE in different stages are deduced. In Section 3, the detailed formulation of ESCOPF model, which consists of the preventive stage, post-contingency stage and corrective stage, is presented. In Section 4, a solving strategy based on Benders decomposition is proposed and the ESCOPF model is effectively solved by an iterative and parallel algorithm. Validities of the proposed model the solving strategy are verified on a 6-bus test system, a modified IEEE RTS-96 system, and a case 2383wp system in Section 5. The main conclusions of this paper are drawn in Section 6.

#### **2. Temperature Variations of Post-Contingency Overload Lines**

The preventive control and corrective control are utilized to guarantee the operation security of power systems in the base case and eliminate long-term power flow violations caused by random failures, respectively. Because corrective measures cannot be implemented without time delay after the contingency, there is a post-contingency stage between the preventive stage and the corrective stage. This time delay comes from two aspects. First, automatic dispatch systems or dispatchers need time to conduct fault locations and determine the fault type. If the automatic reclosing fails, corresponding corrective measures like the unit rescheduling and emergency load shedding should be matched or remade. The time spent in this process is denoted as the response time in this paper. Then, it takes time for conducting those corrective measures such as the unit rescheduling. The specific duration of ramping time depends on adjusted power outputs and ramping rates of units. Therefore, the system post-contingency stage is divided into the response substage and ramping substage while the power flow shows different variations that are presented in Figure 1, and detailed formulations are derived in the following subsections.

**Figure 1.** Two typical power flow and temperature variation curves of the post-contingency overload line *l*: (**a**) the concave curve of the conductor temperature variation; (**b**) the monotonous curve of the conductor temperature variation.

#### *2.1. Preventive Stage*

The preventive stage is the time interval before *t*1. For the discussed line *l*, the power flow *Fl*,0 in this stage is limited by power flow constraints. Since the contingency happens randomly, the conductor temperature of the line *l* in the preventive stage *Tcl*,0 is conservatively assumed to equal its steady-state value. This temperature can be calculated by the steady-state heat balance equation (SHBE) [36].

$$\mathcal{R}\left(T c\_{l,0}\right) I\_{l,0}^2 + q\_5 - q\_\varepsilon (T c\_{l,0}) - q\_\Gamma (T c\_{l,0}) = 0,\tag{1}$$

where *R*(*Tcl*,0) indicates the AC resistance of conductor at the temperature *Tcl*,0 and *Il*,0 denotes the preventive load current which is proportional to the power flow *Fl*,0:

$$I\_{l,0} = I\_{l, \text{Rate}} \frac{F\_{l,0}}{F\_{l, \text{Rate}}} \, \tag{2}$$

where *Il*,*Rate* and *Fl*,*Rate* denote the rated load current and the rated power flow of line *l*, respectively.

In Equation (1), *qs* presents the heat gain rate per unit length from the sun which is regarded as a constant in the short term. *qc* and *qr* are the convected heat loss rate and radiated heat loss rate per unit length of the line *l* which are formulated as functions of *Tcl*,0 in Equations (3) and (4), respectively:

$$q\_c(Tc\_{l,0}) = \left[1.01 + 0.0372 \left(\frac{D\rho\_f V\_W}{\mu\_f}\right)^{0.52}\right] k\_f K\_{\text{angle}} \left(Tc\_{l,0} - Ta\right),\tag{3}$$

$$q\_I(Tc\_{l,0}) = 0.0178D\varepsilon \left[ \left( \frac{Tc\_{l,0} + 273}{100} \right)^4 - \left( \frac{Ta + 273}{100} \right)^4 \right],\tag{4}$$

where *D* indicates the conductor diameter and *VW* denotes the speed of air steam around the line *l*. *ρ<sup>f</sup>* and *μ<sup>f</sup>* are the air density and air dynamic viscosity, respectively. *k <sup>f</sup>* is the thermal conductivity of air. *Kangle* and *Ta* present the wind direction factor and the ambient air temperature, respectively, while *ε* indicates the emissivity. It can be seen from the above two equations that *qc* is the linear function of *Tcl*,0 while the relationship between *qr* and *Tcl*,0 is more complex. Numerical simulations such as the Runge–Kutta method can be utilized to solve the SHBE and obtain *Tcl*,0.

#### *2.2. Response Substage*

After the occurrence of contingency at *t*1, the preventive control still works and the power flow redistributes immediately with the extremely little electric time constant of the system. For the post-contingency overload line, the power transferred on the line *l* increases straight to *Fl*,*kp*(*t*1), which is higher than its rated value. Based on the definition of response time *tresp*, it is expressed as:

$$t\_{resp} = t\_2 - t\_1. \tag{5}$$

The dynamic thermal behavior of the line *l* during this system response substage can be quantified by the transient heat balance equation (THBE) [36].

$$\frac{dT c\_{l,kp}}{dt} = \frac{1}{m \mathbb{C}\_p} \left[ R \left( T c\_{l,kp} \right) I\_{l,kp}^2 + q\_s - q\_c (T c\_{l,kp}) - q\_r (T c\_{l,kp}) \right],\tag{6}$$

where *Tcl*,*kp* indicates the real-time conductor temperature of the line *l* in the post-contingency stage, and *mCp* denotes the total heat capacity of the conductor. In order to obtain the analytical solution of Equation (6), some approximations are needed to be made. First, the resistance-temperature effect is neglected, and *R* is regarded as a constant. Then, the radiated heat loss rate *qr* is locally linearized at the median of the ambient temperature *Ta* and the rated conductor temperature *Tcl*,*Rate*:

$$
\eta\_r(Tc\_{l,kp}) \approx R\_1Tc\_{l,kp} + R\_2(Ta),
\tag{7}
$$

where *R*<sup>1</sup> and *R*<sup>2</sup> are equivalent parameters. THBE can be simplified as a linear differential equation:

$$\frac{dT\mathbf{c}\_{l,kp}}{dt} = K\_1 T\mathbf{c}\_{l,kp} + K\_2 I\_{l,kp}^2 + K\_3 \quad \mathbf{t} \in [\mathbf{t}\_1, \mathbf{t}\_2],\tag{8}$$

where *K*1, *K*<sup>2</sup> and *K*<sup>3</sup> denote parameters depending on the line type and surrounding environment conditions. The variation of load current *Il*,*kp* during this stage is described by a step function:

$$I\_{l,kp}(t) = I\_{l,Rate} \frac{\left(F\_{l,kp}(t\_2) - F\_{l,0}\right)H(t - t\_1) + F\_{l,0}}{F\_{l,Rate}} \quad t \in [t\_1, t\_2],\tag{9}$$

where *H*(*t*) presents the Heaviside step function. The analytical expression of *Tcl*,*kp* are derived from Equations (8) and (9):

$$T\mathbf{c}\_{l,kp}(t) = -\frac{K\_2 I\_{l,kp}^2(t\_1) + K\_3}{K\_1} + \exp\left[\mathbf{K}\_1\left(t - t\_1\right)\right] \left(\frac{K\_2 I\_{l,kp}^2(t\_1) + K\_3}{K\_1} + T\mathbf{c}\_{l,0}\right) \quad t \in [t\_1, t\_2]. \tag{10}$$

It can be observed that the temperature variation of the line *l* in the system response process is an exponential function of the time *t*.

#### *2.3. Ramping Substage*

Corrective measures start to be taken at *t*<sup>2</sup> and the post-contingency power flow *Fl*,*kp* gradually drops below its rated value. At the beginning of ramping substage, all rescheduled units adjust their active power outputs while the power flow of the line *l* decreases fast. After a few mins, some units finish adjustments and the decline becomes slow. Therefore, the curve of actual power flow variation is multiplied piecewise and convex. In order to simplify the problem, it is assumed that *Fl*,*kp* linearly decreases from *Fl*,*kp*(*t*2) to the corrective value *Fl*,*<sup>k</sup>* during the ramping substage. In Figure 1, *t*<sup>3</sup> denotes the time point when the last unit accomplishes its rescheduling plan. The load current *Il*,*kp* in this stage can be approximately formulated as:

$$I\_{l,kp}(t) \approx \frac{I\_{l,Rate}}{F\_{l,Rate}} \left[ \frac{F\_{l,k} - F\_{l,kp}(t\_2)}{t\_{ramp}} (t - t\_2) + F\_{l,kp}(t\_2) \right] \quad t \in [t\_2, t\_3],\tag{11}$$

where the ramping time *tramp* is defined as:

$$t\_{ramp} = t\_3 - t\_2.\tag{12}$$

Equation (11) can be further compacted and represented as:

$$I\_{l,kp}(t) \approx I\_{l,kp}(t\_2) - \Delta I\_{l,kp} \left(t - t\_2\right) \quad t \in [t\_2, t\_3],\tag{13}$$

where *Il*,*kp*(*t*2) and Δ*Il*,*kp* are expressed by the power flow *Fl*,*kp*(*t*2) and *Fl*,*k*:

$$I\_{l,kp}(t\_2) = \frac{I\_{l,Rate}}{F\_{l,Rate}} F\_{l,kp}(t\_2),\tag{14}$$

$$
\Delta I\_{l,kp} = \frac{I\_{l,Rate}}{F\_{l,Rate}} \frac{F\_{l,k} - F\_{l,kp}(t\_2)}{t\_{ramp}}.\tag{15}
$$

*Energies* **2020**, *13*, 48

Because the actual load current of the line *l* at any time from *t*<sup>2</sup> to *t*<sup>3</sup> does not exceed the approximate value provided by Equation (13), the above simplification makes the thermal rating constraint stricter.

The temperature variation can be obtained through Equations (8) and (13):

$$\begin{split} T c\_{l,kp}(t) &= \exp\left[K\_1(t - t\_2)\right] \left(T c\_{l,kp}(t\_2) + K\_4\right) - K\_4 - \frac{\Delta I^2 K\_2 (t - t\_2)^2}{K\_1} - \\ &\quad \frac{2\Delta I K\_2 (t - t\_2) \left(\Delta I\_{l,kp} - I\_{l,kp}(t\_2) K\_1\right)}{K\_1^2} \quad t \in [t\_2, t\_3], \end{split} \tag{16}$$

where

$$K\_{4} = \frac{K\_{2}I\_{l,kp}^{2}(t\_{2})K\_{1}^{2} - 2K\_{2}I\_{l,kp}(t\_{2})\Delta I\_{l,kp}K\_{1} + 2K\_{2}\Delta I\_{l,kp}^{2} + K\_{3}K\_{1}^{2}}{K\_{1}^{3}}.\tag{17}$$

In the initial period of the ramping stage, the real-time conductor temperature *Tcl*,*kp*(*t*) is lower than the steady-state temperature of the real-time power flow *Fl*,*kp*(*t*), hence the temperature continues increasing. There are two different curves of the subsequential temperature variation.

(a) If *Fl*,*kp*(*t*2) is significantly higher than the corrective value *Fl*,*k*, *Fl*,*kp*(*t*) decreases rapidly; there is a stationary point at *thT* when *Tcl*,*kp*(*thT*) equals the steady-state temperature of *Fl*,*kp*(*thT*). The temperature reaches its maximum *Tcmax <sup>l</sup>*,*kp* and then declines with the continuing decrease of the power flow *Fl*,*kp*. This temperature variation is shown in Figure 1a.

(b) If *Fl*,*<sup>k</sup>* is slightly lower than *Fl*,*kp*(*t*2), while *Fl*,*kp*(*t*) declines slowly between *t*<sup>2</sup> and *t*3. The ramping stage may end before the temperature *Tcl*,*kp*(*t*) climbs to the steady-state value of the real-time power flow *Fl*,*kp*(*t*). Therefore, the temperature variation is monotonous in this case, which is presented in Figure 1b.

#### *2.4. Corrective Stage*

After *t*3, corrective measures are finished while the power flow falls below the limit specified by security constraints. The load current of the line *l* is expressed as:

$$I\_{l,k} = I\_{l,Rate} \frac{F\_{l,k}}{F\_{l,Rate}}.\tag{18}$$

The corrective conductor temperature *Tcl*,*<sup>k</sup>* can be derived from Equations (8) and (18):

$$Tc\_{l,k}(t) = -\frac{K\_2 l\_{l,k}^2 + K\_3}{K\_1} + \exp\left[K\_1 \left(t - t\_3\right)\right] \left(\frac{K\_2 l\_{l,k}^2 + K\_3}{K\_1} + Tc\_{l,kp}(t\_3)\right) \quad t \in [t\_{3\prime} + \infty). \tag{19}$$

It can be seen from the above equation that the conduction temperature variation is monotonous. Due to power flow constraints existing in the corrective stage, the security operation of line *l* can be guaranteed as long as *Tcl*,*kp*(*t*3) is lower than the rated temperature.

Based on the electrothermal analysis in the four stages above, the key point of improving the system security is controlling the post-contingency conductor temperature within the permissible scale. The highest temperature could be in the initial period of the ramping substage or at the end of this stage. However, the specific variation type cannot be judged by a succinct algebraic criterion before plotting the curve according to temperature functions. Meanwhile, due to the combination of an exponential function and a quadratic function on the right side of Equation (16), explicit formulations of *thT* and *Tcmac <sup>l</sup>*,*kp* cannot be obtained.

#### **3. Enhanced Security-Constrained OPF**

The enhanced security-constrained OPF model is established on the basis of the PSCOPF model. Thermal rating constraints are integrated in the post-contingency stage to improve the system security and the objective is to minimize the generation cost. Detailed formulations are presented as follows:

$$\min \sum\_{i=1}^{N\_G} \left( a\_i P G\_{i,0}^2 + b\_i P G\_{i,0} + c\_i \right), \tag{20}$$

$$\text{s.t.} \sum\_{i=1}^{N\_G} PG\_{i,0} = \sum\_{j=1}^{N\_D} PL\_{j\prime} \tag{21}$$

$$PG\_{i,min} \le PG\_{i,0} \le PG\_{i,max} \quad i = 1,2,...,N\_{G\_{\prime}} \tag{22}$$

$$-F\_{\text{Rate}} \le F\_0 = T\_0 \left[ PG\_0 - PL \right] \le F\_{\text{Rate}} \tag{23}$$

$$F\_{kp} = T\_k \left[ PG\_0 - PL \right] \quad k = 1, 2, \ldots, N\_{K\_{\prime}} \tag{24}$$

$$\Pr\_{l,kp}^{\text{max}}(F\_{l,0}, F\_{l,kp}, F\_{l,k}, t\_{\text{resp}}, t\_{\text{ramp}}) \le \text{Tc}\_{l,\text{Rate}} \quad l = 1,2,...,N\_{\text{L}} \quad k = 1,2,...,N\_{\text{K}} \tag{25}$$

$$PG\_{i,min} \le PG\_{i,k} \le PG\_{i,max} \quad i = 1,2,...,N\_G \quad k = 1,2,...,N\_K \tag{26}$$

$$PG\_{i,0} - \Lambda PG\_i \le PG\_{i,k} \le PG\_{i,0} + \Lambda PG\_i \quad i = 1, 2, \dots, N\_G \quad k = 1, 2, \dots, N\_{\mathcal{K}} \tag{27}$$

$$1 - F\_{\text{Rate}} \le F\_{\text{k}} = T\_{\text{k}} \left[ PG\_{\text{k}} - PL \right] \le F\_{\text{Rate}} \quad k = 1, 2, \dots, N\_{\text{K}}.\tag{28}$$

In the objective function (20), *ai* and *bi* denote the second order and the first order cost coefficients of unit *i*, respectively. *ci* indicates the fixed cost. *PGi*,0 is the active power output of the unit *i* determined by the preventive control. *NG* is the number of units. Because probabilities of contingencies are quite low and the rescheduling amount of power output is limited by unit ramping capacities, the cost of corrective control can be neglected.

Constraints (21)–(23) are utilized to guarantee the system operation security in the preventive stage. Equation (21) implies the active power balance where *PLj* is the load demand at the bus *j* and *NL* represents the number of load buses. Constraints (22) are power output limits where *PGi*,*min* and *PGi*,*max* denote the minimum and maximum power output of the unit *i*, respectively. In power flow constraints (23), *PG***<sup>0</sup>** and *PL* are the preventive power output vector and the load demand vector, respectively. *FRate* and *F***<sup>0</sup>** indicate the rated power flow vector and the preventive power flow vector, respectively. *T***<sup>0</sup>** is the shift matrix in the base case without any failures. In the post-contingency stage, due to the change of grid topology caused by device outages, the power flow distribution is needed to be recalculated for each contingency. Constraint (24) indicates the power flow equation where *Fkp* and *Tk* are post-contingency power flow vector and the shift matrix in the contingency *k*, respectively. Conductor temperatures are limited by the constraint (25). According to the analysis of temperature variation in Section 2, the highest temperature *Tcmax <sup>l</sup>*,*kp* in the post-contingency stage *k* depends on *Fl*,0, *Fl*,*kp*, *Fl*,*k*, *tresp* and *tramp*, and their functional relationships are determined by Equations (5)–(17).

Corrective control constraints consist of inequations (26)–(28). Constraints (26) and (27) indicate unit output limits in the corrective stage, where *PGi*,*<sup>k</sup>* is the rescheduled power output of the unit *i* in the contingency *k* and Δ*PGi* is the adjustable power output determined by its reserve capacity. In the power flow constraint (28), *Fk* denotes the corrective power flow vector and *PGk* indicates the rescheduled power output vector in the contingency *k*. Line overloads must be completely eliminated by corrective measures and the power system enters a new safe operation status.

The established ESCOPF model is a multi-stage dispatch problem with the quadratic objective function (20) and nonlinear constraints (25). Due to the complex form of the temperature variation function in Equation (16), *Tcmax <sup>l</sup>*,*kp* is presented as an implicit function about *Fl*,0, *Fl*,*kp*, *Fl*,*k*, *tresp* and *tramp*. This implies that conventional optimization methods could not be directly utilized to solve the established model, and some pretreatment measures need to be made. Moreover, different unit rescheduling plans are made for each contingency scenario. Therefore, the computational burden increases as the number of contingencies grows. It is necessary to control the solving time when applying the proposed model in large-scale power systems

#### **4. Solving Strategy Based on Benders Decomposition**

Benders decomposition is one of the commonly used decomposition techniques in power systems. J. F. Benders first introduced this decomposition algorithm for solving large-scale mixed-integer programming (MIP) problems [37]. At present, it has been successfully applied to handle various optimization problems including the unit commitment [38], economic dispatch [39], and transmission planning [40]. In applying Benders decomposition, the optimal solution may require iterations between the master problem and several subproblems. The lower bound and upper bound solution are provided by the master problem and feasible subproblems, respectively. If any of the subproblems is infeasible, a Benders cut will be introduced to the master problem. When the upper bound and the lower bound are sufficiently close, the iteration stops. In this section, a solving strategy based on Benders decomposition is presented to solve the established model with nonlinear thermal rating constraints.

The ESCOPF model is first decomposed into a master problem optimizing the preventive control and two subproblems for the corrective control feasibility check and the line thermal rating check, respectively. In the second subproblem, conductor temperatures are calculated based on the determined power flow distribution in the preventive, post-contingency and corrective stage. Meanwhile, the independence between each contingency is taken into account and both subproblems can be processed in parallel. For those contingencies which cannot pass the corrective feasibility check or the line thermal rating check, corresponding Benders cuts are generated and returned into the master problem. The master problem and subproblems are sequentially solved until both checks are satisfied for all contingency scenarios. The specific solving flow is shown in Figure 2 and detailed formulations for the master problem and two subproblems are presented as follows.

**Figure 2.** Solving flow based on Benders decomposition.

#### *4.1. Master Problem (Preventive Control Optimization)*

The preventive control problem is optimized separately in the master problem which consists of the objective function (20), constraints (21)–(23). The preventive power flow is obtained and the power flow redistribution in each contingency during the response stage can be calculated by Equation (24). Corresponding temperature variations from *t*<sup>0</sup> to *t*<sup>2</sup> are quantified by Equations (1) and (10). In each iteration, Benders cuts (29) and (30) are generated for those check violating contingencies and added into the master problem as constraints:

$$\text{Corrective control feasibility check Benders cuts : Constraint} \left( 37 \right), \tag{29}$$

Thermal rating check Benders cuts : Constraints (38). (30)

The feasible region of the master problem is continuously shrunk by returned Benders cuts and the most economical operating point is searched again in the new feasible region. According to following generating approaches for Benders cuts, above constraints (29) and (30) are linear. Therefore, the master problem in each iteration is a linearly constrained quadratic optimization problem which can be directly solved by current optimizers such as Cplex, Gurobi and so on.

#### *4.2. Subproblem 1 (Corrective Control Feasibility Check)*

Due to the existence of correlated constraints (27) between the preventive control and corrective control, the feasibility of corrective control needs to be checked for the operating point determined by the master problem in each iteration. In the credible contingency *k*, slack variables *si*,*<sup>k</sup>* and *ri*,*<sup>k</sup>* are introduced for power output rescheduling constraints. Detailed formulations of the corrective control feasibility check subproblem are presented as follows:

$$\min \sum\_{i=1}^{N\_G} (s\_{i,k} + r\_{i,k}), \tag{31}$$

$$\text{1s.t.}\quad \text{Constrraints} \left(26\right) \text{and} \left(28\right), \tag{32}$$

$$PG\_{i,k} - PG\_{i,0}^\* - s\_{i,k} \le \Delta PG\_i \quad i = 1, 2, \dots, N\_{G\_{\prime\prime}} \tag{33}$$

$$PG\_{i,k} - PG\_{i,0}^\* + r\_{i,k} \ge -\Delta PG\_i \quad i = 1, 2, \dots, N\_{\mathbb{G}^\star} \tag{34}$$

$$s\_{i,k} \ge 0 \quad i = 1, 2, \dots, N\_{G'} \tag{35}$$

$$r\_{i,k} \ge 0 \quad i = 1, 2, \dots, N\_{\mathbb{G}}.\tag{36}$$

In the constraints (33) and (34), the determined preventive power output *PG*∗ *<sup>i</sup>*,0 is a constant provided by results of the master problem in the current iteration. After solving the subproblems (31)–(36), the corrective power flow distribution is obtained and temperature variations after *t*<sup>2</sup> can be quantified by Equations (16) and (19).

Once the optimized result of the objective function (31) equals 0, overflows in the contingency *k* can be eliminated by a feasible unit rescheduling. Otherwise, corresponding Benders cuts are needed to be generated and returned to the master problem. The Benders cut of the contingency *k* which cannot be corrected is deduced as:

$$\sum\_{i=1}^{N\_G} \left( s\_{i,k} + r\_{i,k} \right) + \sum\_{i=1}^{N\_G} \left[ \left( -\lambda\_{i,k} + \gamma\_{i,k} \right) \left( PG\_{i,0} - PG\_{i,0}^\* \right) \right] \le 0,\tag{37}$$

where *λi*,*<sup>k</sup>* and *γi*,*<sup>k</sup>* are the Lagrange multipliers of the constraint (33) and (34), respectively. The Benders cut (37) is a linear constraint about *PGi*,0 where *si*,*k*, *ri*,*k*, *λi*,*k*, *γi*,*<sup>k</sup>* and *PG*<sup>∗</sup> *<sup>i</sup>*,0 are constants. Only if the optimized result of this subproblem equals 0 for all contingencies, the operating point set by the master problem can pass the feasibility check for the corrective control.

#### *4.3. Subproblem 2 (Line Thermal Rating Check)*

Temperature variations for four sequential stages are calculated based on the power flow variation determined by the master problem and the corrective control feasibility check subproblem. In each contingency, only if the highest temperature of each line conductor is below its rated value, the thermal rating check can be satisfied. Otherwise, the Benders cut is needed to be derived and added into the master problem:

$$\sum\_{l \in L\_k} \left( T\_{l,kp}^{\text{max}} - T c\_{l,Ratter} \right) + \sum\_{i=1}^{N\_G} \left( \sum\_{l \in L\_k} \left. \frac{\partial T c\_{l,kp}^{\text{max}}}{\partial P G\_{i,0}} \right|\_{P G\_{i,0} = P G\_{i,0}^\*} \right) \left( P G\_{i,0} - P G\_{i,0}^\* \right) \le 0,\tag{38}$$

where *Lk* denotes the set of transmission lines which violate thermal rating constraints. The key point of obtaining the thermal rating Benders cut is calculating the partial derivative of *Tcmax <sup>l</sup>*,*kp* with respect to *PGi*,0. In this paper, an equivalent time method is proposed to build an explicit relationship between the highest conductor temperature and the preventive power output.

The basic idea is using an extended response process to simulate the original two-stage process of the post-contingency temperature variation. The equivalent principle is that those two processes have the same highest temperature and the equivalent time *t equal <sup>l</sup>*,*kp* of the line *l* in the contingency *k* can be calculated from the following equation:

$$\mathrm{Tc}\_{l,kp}^{\mathrm{max}} = -\frac{K\_2 I\_{l,kp}^2(t\_1) + K\_3}{K\_1} + \exp\left(\mathrm{K}\_1 \mathrm{t}\_{l,kp}^{\mathrm{quad}}\right) \left(\frac{K\_2 I\_{l,kp}^2(t\_1) + K\_3}{K\_1} + \mathrm{Tc}\_{l,0}\right). \tag{39}$$

Due to the monotonous variation of the conductor temperature in the extended response process, the partial derivative of *Tcmax <sup>l</sup>*,*kp* with respect to *PGi*,0 can be derived by the chain rule:

$$\frac{\partial \mathcal{I} \mathbf{C}\_{l,kp}^{\text{max}}}{\partial \mathcal{I} \mathbf{G}\_{l,0}} = \left\{-1 + \exp\left(\mathcal{K}\_{l} \mathbf{f}\_{l,kp}^{\text{equal}}\right)\right\} \frac{\mathcal{K}\_2}{K\_1} 2\mathcal{I}\_{l,kp}^\* \frac{\partial \mathcal{I}\_{l,kp}(t\_1)}{\partial \mathcal{I} \mathbf{G}\_{l,0}} + \exp\left(\mathcal{K}\_l \mathbf{f}\_{l,kp}^{\text{equal}}\right) \frac{\partial \mathcal{I} \mathbf{c}\_{l,0}}{\partial \mathcal{I} \mathbf{G}\_{l,0}}\tag{40}$$

where *I*∗ *<sup>l</sup>*,*kp* indicates the load current of the line *l* at *t*<sup>1</sup> in the contingency *k* which can be obtained from results of the master problem. The partial derivative of *Il*,*kp* with respect to *PGi*,0 can be calculated from the power flow Equation (24):

$$\frac{\partial I\_{l,kp}}{\partial P G\_{l,0}} = \frac{dI\_{l,kp}}{dF\_{l,kp}} \frac{\partial F\_{l,kp}}{\partial P G\_{l,0}} = \frac{I\_{l,R, \text{Rate}}}{F\_{l,R \text{Rate}}} T\_k(l,i). \tag{41}$$

Similarly, based on the power flow equation in the preventive stage (23), the partial derivative of *Tcl*,0 with respect to *PGi*,0 can be deduced as:

$$\frac{\partial T\mathbf{c}\_{l,0}}{\partial P\mathbf{G}\_{i,0}} = \frac{d T\mathbf{c}\_{l,0}}{d I\_{l,0}} \frac{d I\_{l,0}}{d F\_{l,0}} \frac{\partial F\_{l,0}}{\partial P\mathbf{G}\_{i,0}} = \frac{d T\mathbf{c}\_{l,0}}{d I\_{l,0}} \frac{I\_{l,Rate}}{F\_{l,Rate}} T\_0(l,i). \tag{42}$$

Using a quadratic function to fit the SHBE, the above partial derivation can be further derived as:

$$T\mathbf{c}\_{l,0} = M\_1 I\_{l,0}^2 + M\_2 I\_{l,0} + M\_3 \,\tag{43}$$

$$\frac{dT c\_{l,0}}{dI\_{l,0}} = 2M\_1 I\_{l,0}^\* + M\_2. \tag{44}$$

In Equation (44), *I*∗ *<sup>l</sup>*,0 denotes the determined preventive load current of the line *l*. It can be seen from the Equations (41)–(44) that the partial derivative of *Tcmax <sup>l</sup>*,*kp* with respect to *PGi*,0 is a constant for the determined operating point and the unit rescheduling plan. This means the return Benders cut (38) is also a linear constraint.

The pseudocode of the solving strategy based on Benders decomposition is shown in Algorithm 1. It may help to reproduce and implement the approach proposed in this paper. The *NFC* and *NTC* indicate the number of contingenies which cannot pass the corrective control feasibility check and the line thermal rating check.

#### **Algorithm 1** The solving strategy based on Benders decomposition

**Input:** Power system parameters and environment parameters. Durations of the response time and

```
ramping time. The rated conductor temperature. The contingency set.
 1: while NFC = 0 and NTC = 0 do
 2: solve the master problem consisting of Equations (20), (21)–(23) and (29), (30);
 3: obtain unit active output PG∗
                               0 and power flow distribution F0;
 4: NFC = 0;
 5: for each k ∈ [1, NK] do
 6: solve the subproblem 1 consisting of Equations (31)–(36);
 7: if ∑ si,k + ri,k = 0 then
 8: generate the Benders cut by Equation (37) and return it into th master problem;
 9: obtain the power flow distribution Fk; NFC = NFC + 1;
10: end if
11: end for
12: NTC = 0;
13: for each k ∈ [1, NK] do
14: for each line l operating in the contingency k do
15: calculate the Tcmax
                         l,kp and t
                                 equal
                                 l,kp by Equations (1)–(19) and (39), repsectively;
16: if Tcmax
               l,kp > Tcl,Rate then
17: add line l into the set Lk
18: end if
19: end for
20: if Lk = φ then
21: NTC = NTC + 1; sumdTc = 0;
22: for each l ∈ Lk do
23: sumdTc = sumdTc + (Tcmax
                                   l,kp − Tcl,Rate);
24: end for
25: for each i ∈ [1, NG] do
26: sumpdi = 0;
27: for each l ∈ Lk do
28: calculate the partial derivative
                                            ∂Tcmax
                                               l,kp
                                             ∂PGi,0 by Equations (40)–(44);
29: sumpdi = sumpdi +
                                   ∂Tcmax
                                     l,kp
                                   ∂PGi,0 ;
30: end for
31: end for
32: return sumdTc + NG
                           ∑
                          i=1
                             [sumpdi(PGi,0 − PG∗
                                              i,0))] ≤ 0 into the master problem;
33: end if
34: end for
35: end while
Output: The active power output of each unit in the preventive stage and its adjusted amount in the
```
corrective stage. The power flow distribution and variations of line conductors.

#### **5. Case Studies**

The model and solving strategy proposed in this paper are validated on a 6-bus test system, a modified IEEE RTS-96 system, and a case 2383wp test system. It assumed that all transmission lines in those three test systems use the 400 mm<sup>2</sup> Drake 26/7 ACSR conductor (Nexans S.A., Paris, France) whose parameters (resistance *R*, diameter *D*, heat capacity *mCp* and emissivity *ε*) are given by Reference [30] and shown in Table 1.


**Table 1.** The parameters of transmission lines.

Another assumption is that lines operate in the same environment condition, although it is difficult to be satisfied for large-scale power systems. The related environment parameters can be divided into two kinds. The air density *ρ<sup>f</sup>* , dynamic viscosity *μ<sup>f</sup>* and the thermal conductivity of air *k <sup>f</sup>* may not have obvious variations and could be regarded as constants. The other parameters including the ambient air temperature *Ta*, speed of air steam *VW*, wind direction factor *Kangle* and the power gain rate from the sun *qs* can be provided by online micro-meteorology monitoring systems. With the development of smart grids, those systems have a wide range of applications [41]. The more accurate measured parameters are obtained, the more effective control for the line thermal rating can be achieved. The data collected by distributed monitoring systems are transmitted through the high-speed optical networks or wireless networks like ZigBee, GPRS and 4G to energy management systems (EMS). Through data preprocessing, those data can be used to solve the electro-thermal coupling OPF model. Values of two kinds of parameters are presented in Table 2.

**Table 2.** The environment parameters.


In addition, for protection systems, the highest threshold value of line load current should be appropriately raised and the time-delay of triggering for overload lines needs to be increased. It will let protections not take action immediately and allow transmission lines to operate under the overload condition for a short period of time. In addition, then the ESCOPF model proposed in this paper can work in practical applications.

Under the condition of given parameters, the time constant of line thermal process is about 14 min which has been verified in Reference [30]. Meanwhile, the rated load current of lines is set at 992 A while the corresponding steady-state temperature is 100.0 ◦C. The widely used "N−1" criterion is adopted for the former two test systems to filter line contingencies. This means, in each contingency, there is just a single transmission line outages. All simulations are performed on a personal computer with 4 Intel (R) Core (TM) i5-6200U CPU (2.3 GHz) (Lenovo, Beijing, China) and 8 GB memory. The programs are implemented using a Matlab (R2016a, The MathWorks, Natick, MA, USA) environment and underlying optimization problems are solved by Cplex.

#### *5.1. 6-Bus Test System*

The 6-bus test system consists of three units, three load bus, and 11 transmission lines. Its wiring diagram is presented in Figure 3. Detailed parameters of units, load buses and lines are given by Tables A1–A3 (Appendix A), respectively. There are 11 "N−1" contingencies where the line *k* failures in the contingency *k*. The response time and the ramping time of this system are 5 min and 7 min, respectively, and the rated conductor temperature is conservatively set at 100.00 ◦C. The mean program executing time of 10 tests with the consistent convergence solution is 10.07 s. Simulation results in each iteration are shown in Table 3.

In the initial iteration, the preventive control is optimized separately with the lowest generation cost \$861.92. Four contingencies have post-contingency overflows which cannot be eliminated by the feasible corrective control. The thermal rating check is violated in three contingencies, and the highest temperature reaches 136.64 ◦C. Unit active power outputs are adjusted by returned Benders cuts as the iteration number grows. Combining with generation cost coefficients in Table A1, it can be found that the output of the most expensive unit G2 increases, which results in the growth of *TGC*. It is worth noting that, after iteration 3, *NFC* increases from 0 to 1, while *NTC* decreases from 2 to 1. This implies that the corrective control feasibility check and the thermal rating check may be conflicting in some situations, and this conflict is solved by further shrinking the feasibility region of the operating point. After five iterations, both *NFC* and *NTC* drop to 0 and the highest conductor temperature is strictly controlled within the permissible scale 100.00 ◦C. This means that both short-term and long-term safe operation of this 6-bus system can be guaranteed. The final generation cost rises to \$917.51, which is 6.5% higher than its initial results.

**Figure 3.** The wiring diagram of 6-bus test system.


**Table 3.** The simulation results in each iteration for a 6-bus system.

*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>c</sup>* The total generation cost.

In addition, temperature variations of lines that have the highest conductor temperature in each contingency are shown in Figure 4. Displayed durations of the preventive stage and the corrective stage are both 5 min. The post-contingency stage starts at the 5th min and ends at the 17th min. It is further confirmed that conductor temperatures in all contingencies are below the allowable value 100.00 ◦C. The highest temperature occurs on line 1, which connects two units G1 and G2 in the "N−1" contingency with the failure of line 2. The post-contingency load rate of line 1 reaches 1.35 at the 5th min and drops to 1.00 from the 10th min to the 17th min. The temperature rapidly rises to the peak 100.00 ◦C at around the 15th min and then slowly declines. This result verifies that the power flow constraint and the thermal rating constraint are not equivalent, and the former constraint is over conservative from the short-term perspective.

The conductor temperature varies monotonously under the condition that the corrective power flow is higher or slightly lower than it in the post-contingency stage. This corresponds to cases of the highest temperature line in the contingencies 1, 3, 5, 6, and 7. For those lines with highest temperature in the contingencies 2, 4, 8, 9, 10, and 11, the post-contingency power flow is significantly higher than the corrective value while concave temperature curves are presented. In those cases, the maximum absolute difference between load rates in the post-contingency and the corrective stage reaches 0.35, which is obviously higher than the value 0.14 in former cases. This result validates the analysis about two typical temperature variations in Section 2.

**Figure 4.** The temperature variations of the line with the highest temperature in each contingency.

It can be found from Figure 4 that the highest temperature occurs five times on line 9 in 11 credible contingencies. Line 9 directly connects the unit G3 and the load bus L3 while the load rate in the preventive stage is 0.90. It is markedly higher than the second highest load rate 0.75 of line 3, which becomes the highest temperature line for the remaining two contingencies. From this perspective, line 9 can be regarded as the key transmission line in the 6-bus system while the security and economy of the system could be improved by appropriately raising its transmission capacity.

#### *5.2. IEEE RTS-96 System*

The modified RTS-96 system consists of three interconnected RTS-79 systems and has 96 units, 73 buses, and 120 branches in total. Its topology is presented in Figure 5. The unit data and load data can be found in Reference [42] while the modified line data are provided by Reference [15]. Because line 52 in area 2 and line 90 in the area 3 are single-circuit lines connecting to the bus 207 and bus 307, respectively, their outages will cause a power imbalance in the post-contingency stage. Other measures

such as more frequent inspections and maintenance could be used to improve their secure operations. Therefore, the final number of "N−1" line contingencies is 118.

**Figure 5.** The wiring diagram of the IEEE RTS-96 test system.

In the base case, the response time and ramping time of this system are also set at 5 min and 7 min, respectively, while the rated conductor temperature is 100.00 ◦C. The mean value of the program executing time for 10 simulations with the same convergence solution equals 53.87 s and detailed simulation results in each iteration are shown in Table 4.


**Table 4.** The simulation results in each iteration for the IEEE RTS-96 system.

*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>c</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>d</sup>* The total generation cost.

There are 10 contingencies which cannot satisfy the corrective control feasibility check in the initial solution. Post-contingency overflows occur in all contingencies, and the highest load rate reaches 1.78. From the perspective of thermal rating, 19 contingencies are not secure. The highest temperature is 136.04 ◦C, and the lowest total generation cost is \$5058.8. After the first iteration, *NFC*, *NPC*, and *NTC* decrease significantly to 2, 53, and 8, respectively. In results after iteration 3, *TGC* approaches its final optimized value. This is because the same type of units have the same generation cost coefficients. Actually, the generation plan is still adjusted, which can be observed from variations of *NFC*, *LRmax*, *NTC*, and *Tcmax*. After seven iterations, long-term overflows can be totally removed by the feasible corrective control while the highest temperature equals the rated value 100 ◦C, which confirm that the algorithm has converged. It is worth noting that there are still 51 contingencies having post-contingency overload lines. However, according to the analysis in Section 2, those short-term power flow violations will not affect the operation security of transmission lines.

In order to further verify the effective control for conductor temperatures through the ESCOPF model, the highest conductor temperature of each transmission line in all credible contingencies is presented in Figure 6.

**Figure 6.** The highest conductor temperature of each line in all credible contingencies.

It can be seen that most of temperatures are in the interval between 50 ◦C and 90 ◦C. The highest post-contingency temperature 100 ◦C occurs on line 10 and line 11, while temperatures of line 51, line 53, line 89, and line 91 are also very close to this limit. Based on the system diagram, it can be found that those lines are at similar topological locations in corresponding areas. For instance, line 10, line 51, and line 89 connect the bus 106 and bus 110, bus 206 and bus 210, bus 306 and bus 310, respectively. When focusing on the area 1, the 136 MW load demand is located at the bus 106. Once line 5 connecting bus 102 and bus 106 fails, 136 MW of power is needed to be transferred by line 10 alone. Meanwhile, the load rate of line 10 in the preventive stage 0.82 is relatively higher. The system operation security can be improved by controlling conductor temperatures of those key transmission lines in the post-contingency stage.

The performance of the proposed ESCOPF model is further compared with the conventional PSCOPF model and CSCOPF model from the perspective of security and economy, while detailed results are listed in Table 5.


**Table 5.** The comparison of results of PSCOPF, ESCOPF, and CSCOPF models.

*<sup>a</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>b</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>c</sup>* The total generation cost.

It can be seen that the PSCOPF is the most conservative model with the highest generation cost \$5145.4. For the operating point obtained through this model, power flow constraints are satisfied for all contingencies. Since there is no post-contingency stage, the values of *NTC* and *Tcmax* do not exist. The CSCOPF model tries to eliminate overflows by the corrective control while effects of the time delay are ignored. Both power flow constraints and thermal rating constraints are violated in the post-contingency stage. The highest load rate and conductor temperature are 1.49 and 115.21 ◦C, respectively. Total generation cost \$5063.6 is the lowest. In results of the proposed ESCOPF model, *NPC* declines from 88 to 51 and *LRmax* is limited to 1.34. All credible contingencies can pass the thermal rating check. From the viewpoint of line dynamic thermal behavior, the ESCOPF model and the PSCOPF model have the equivalent security. However, the former total generation cost \$5072.5 is just slightly higher than the result of CSCOPF model.

In addition, influences of the rated conductor temperature on the highest line load rate in all contingency scenarios and the total generation cost of the system are studied and results are shown in Figure 7.

**Figure 7.** The highest load rate and total generation cost under various rated temperatures: (**a**) the highest load rate variation; (**b**) the total generation cost variation.

The highest load rate in the post-contingency stage increases while the total generation cost decreases as the rated conductor temperature grows. It can be observed that both rates of change decline. When the rated temperature is set at 88.53 ◦C, there is no contingency having overload lines, and the total generation cost equals the result of the PSCOPF model. Actually, the corrective control is no more needed while the ESCOPF model degenerates into a PSCOPF model. When the preset rated temperature is higher than 115.21 ◦C which is over the highest temperature occurring in results of the CSCOPF model, thermal rating constraints in the post-contingency stage become redundant. This means that the established ESCOPF model is converted into a CSCOPF model and their final optimized generation costs \$5063.6 are consistent.

Moreover, Figure 8 presents the surface of total generation costs under various combinations of the response time and the ramping time. The generation cost remains unchanged while the response time is below 3 min and the ramping time is within 5 min. Due to the thermal inertia, lines can operate under extremely high load rates and do not violate thermal rating limits for a short duration. The operating point of the system is not affected by post-contingency thermal rating constraints. As both growths of the response time and the ramping time, the generation cost continues increasing. It can be observed that the variation of the response time has a greater influence. This is because the post-contingency power flow keeps up a high level in this duration. While the response time is in the interval of 5–8 min and the ramping time is within the interval of 6–9 min, the increase of the total generation cost is the most significant. Beyond this time interval, the growth becomes slow. The above analysis implies that shortening the system response time and implementing units with the quick ramping capacities can improve the security of power systems and decline its operating cost.

**Figure 8.** The total generation cost under various combinations of the response time and ramping time.

#### *5.3. Case 2383wp Test System*

In order to verify the validities of the proposed ESCOPF model and solving strategy on large-scale power systems, the case 2383wp test system that consists of 2383 buses, 2896 branches, and 327 units is adopted. Its detailed system parameters can be found in the file provided by Matpower (6.0, Cornell University, Ithaca, NY, USA). In addition, the "N−5" contingencies are randomly generated for simulations. It means that there are five failure transmission lines in each contingency. The rated conductor temperature is set at 100.00 ◦C for all lines. Simulations are conducted for 10, 20, 50 and 100 "N−5" contingencies and the executing time is 367.83 s, 547.76 s, 1817.73 s, and 4540.95 s, respectively. The detailed results are shown in Table 6–9.


**Table 6.** The simulation results for case 2383wp test system under 10 "N−5" contingencies.

*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>c</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>d</sup>* The total generation cost.

**Table 7.** The simulation results for case 2383wp test system under 20 "N−5" contingencies.


*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>c</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>d</sup>* The total generation cost.

**Table 8.** The simulation results for case 2383wp test system under 50 "N−5" contingencies.


*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>c</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>d</sup>* The total generation cost.

It can be found that the algorithm can converge quickly and steadily for the four simulations above within 8 iterations. The executing time increases with the growth of the number of contingencies. Those results are obtained on a person computer with series computation. According to the analysis in Section 4, both subproblems can be processed in parallel due to the independence between each

contingency. Therefore, when the model and algorithm are applied on more large-scale contingency sets, the parallel computing architecture can be used to effectively decline the executing time.


**Table 9.** The simulation results for case 2383wp test system under 100 "N−5" contingencies.

*<sup>a</sup>* The number of contingencies which cannot pass the corrective control feasibility check, *<sup>b</sup>* The number of contingencies which violate power flow limits in the post-contingency stage, *<sup>c</sup>* The number of contingencies which cannot pass the thermal rating check, *<sup>d</sup>* The total generation cost.

In initial solutions, maximum load rates reach around 8.00 and corresponding highest conductor temperatures are over 1200 ◦C. Both indicators are gradually controlled and decrease significantly as the iterations progress. In final convergence solutions, maximum load rates decline to around 1.6 while the highest temperatures are below the rated value 100.0 ◦C. For optimized operating points, all contingencies can pass the corrective control feasibility check and line thermal rating check. It also needs to be noted that the total generation cost increases from \$27,027 to \$33,593 with the growth of the number of "N−5" contingencies. This is because the feasible region of the power system is shrunk by adding more credible contingencies.

#### **6. Conclusions**

In this paper, an ESCOPF model where thermal rating constraints are integrated to limit post-contingency conductor temperatures is established. After the contingency occurrence, the system first experiences a response stage and a ramping stage in time sequence. Due to the thermal inertia, lines can transfer the power flow that is much higher than the rated value while conductor temperatures are still within the safe scale. The corrective control is implemented with a time delay and long-term overflows are eliminated.

A solving strategy based on Benders decomposition is proposed to deal with the ESCOPF model. The original dispatch problem is divided into a master problem and two subproblems. The system operating point is determined in the master problem while the corrective control feasibility check and the line thermal rating check are separately conducted in two subproblems. The partial derivative of the highest temperature with respect to unit power output in the thermal rating Benders cut is calculated by a presented equivalent time method.

Simulation results on a 6-bus test system, a modified IEEE RTS-96 system, and a case 2383wp test system demonstrate that proposed approaches can offer the following advantages:

(a) The solving strategy based on Benders decomposition steadily converges within eight iterations for three test systems. The finally obtained operating point can pass the feasibility check and thermal rating check and has the lowest generation cost.

(b) The explicit relationship between the highest temperature and unit power output can be built by the proposed equivalent time method. Post-contingency temperatures of line conductors are strictly controlled below the prescribed limit.

(c) The proposed ESCOPF model can find a better balance between the security and economy compared with the conventional SCOPF model. As rated temperatures decline, the model gradually degenerates into the PSCOPF model. Conversely, the obtained operating point approaches the result of the CSCOPF model.

(d) The generation cost increases rapidly as the duration of the response stage rises. Shortening the system duration in the response stage can significantly extend the safe operation region of systems and decrease the total generation cost.

The proposed ESCOPF model is based on the DC power flow which will cause the underestimation of actual conductor temperatures of transmission lines and a more accurate model is needed to be established. In addition, just like the conventional PSCOPF model and CSCOPF model, the ESCOPF model proposed in this paper is a tool that can be used to make the day-ahead, real-time generation plan, and so on. Therefore, a lot of work such as applying this model into the electricity market operation and power systems with the larger-scale renewable energy integration need to be done in future research.

**Author Contributions:** Conceptualization, C.G.; Methodology, X.L.; Validation, W.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work conducted by Xiansi Lou, Wei Chen, and Chuangxin Guo was supported in part by the National Key R&D Program of China under Grant 2017YFB0902600 and in part of the State Grid Corporation of China Project under Grant SGJS0000DKJS1700840.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**


**Table A1.** The parameters of units.


**Table A2.** The parameters of load buses.

**Number "From" Bus Number "To" Bus Number Resistence (p.u.) Rated Power Flow (MVA)** 1 1 2 0.20 50 2 1 4 0.20 70 3 1 5 0.30 55

4 2 3 0.25 55 5 2 4 0.10 80 6 2 5 0.30 40 7 2 6 0.20 70 8 3 5 0.26 50 9 3 6 0.10 80 10 4 5 0.40 40 11 5 6 0.30 40


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
