Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Two-layer Data-driven Robust Scheduling for Industrial Heat Loads  PDF

  • Chuanshen Wu (Member, IEEE)
  • Yue Zhou (Member, IEEE)
  • Jianzhong Wu (Fellow, IEEE)
School of Engineering, Cardiff University, Cardiff CF24 3AA, U.K.

Updated:2025-01-22

DOI:10.35833/MPCE.2024.000105

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

This paper establishes a two-layer data-driven robust scheduling method to deal with the significant computational complexity and uncertainties in scheduling industrial heat loads. First, a two-layer deterministic scheduling model is proposed to address the computational burden of utilizing flexibility from a large number of bitumen tanks (BTs). The key feature of this model is the capability to reduce the number of control variables through analyzing and modeling the clustered temperature transfer of BTs. Second, to tackle the uncertainties in the scheduling problem, historical data regarding BTs are collected and analyzed, and a data-driven piecewise linear Kernel-based support vector clustering technique is employed to construct the uncertainty set with convex boundaries and adjustable conservatism, based on which robust optimization can be conducted. The case results indicate that the proposed method enables the utilization of flexibility in BTs, improving the level of onsite photovoltaic consumption and reducing the aggregated load fluctuation.

I. Introduction

IN the past decades, dispatchable loads have attracted vital attention due to their substantial potential for balancing power systems through management. In industrial sectors, industrial heat loads (IHLs) possess inherent operational flexibility due to their thermal storage capacity, enabling them to participate in demand response in power systems [

1]. During the production processes of IHLs, a considerable amount of energy is required to maintain them within an appropriate temperature range for use at any time. Therefore, the efficient management suggests a significant opportunity for power system balancing [2].

In this paper, an industrial site with bitumen tanks (BTs) is studied as an example of IHLs. The global bitumen market was valued at around 105 billion dollars in 2022, underscoring its substantial presence in the construction and infrastructure sectors worldwide [

3]. In the UK, there are nearly 300 bitumen plants manufacturing over 20 million t of bitumen annually [4]. To ensure the fluidity of bitumen for use at any time, BTs are needed to maintain the bitumen temperature within the required range. At present, there are mainly three heating methods for BTs, namely fuel heating, thermal oil heating, and electric heating [5]. Electric heating utilizes an electric resistance assembly to heat the bitumen in the tank, which is widely applied, especially in applications where precise temperature control of bitumen is essential. According to a survey conducted by the MBA Asphalt Plants [6], one of the largest companies in the bitumen industry, electric BTs account for approximately 15% of the total market share of BTs worldwide, with around 65% of them being used in European countries. Moreover, considering that electric BTs are environmentally friendly due to the absence of combustion of fossil fuels, the proportion of electric BTs is anticipated to continue rising in the low-carbon-oriented future.

Electric BTs are typically equipped with electric heating devices to maintain the bitumen within a certain temperature range, ensuring its fluidity for industrial applications in road construction, building infrastructure, petroleum, chemical, and other various fields [

7]. By properly controlling the on/off status of electrical switches, BTs act as thermal energy storage, providing the flexibility to respond to demand for the power grid. In recent years, there have been studies exploring the possibility of leveraging BTs to provide demand response services. For instance, researchers have developed control algorithms to alter the power consumption of BTs for real-time balancing between supply and demand in electric power systems [8]. Similar methods have been applied to BTs for participating in enhanced frequency response within the Great Britain power system [9]. Existing studies have verified the significant potential of BTs to provide demand response and ancillary services to the power grid.

Nowadays, for tapping the demand response potential of IHLs, there are mainly two research challenges. First, the scheduling of IHLs presents a significant computational challenge due to the strong nonlinearity and large number of integer variables involved [

10]. This is a non-deterministic polynomial-time hard (NP-hard) problem with significant calculational complexity, and thus challenging. Second, the heat transfer rate of IHLs is nonlinear over time and easily affected by uncertainties, such as equipment parameters and environmental conditions [11]. Moreover, both model and data uncertainties propagate over time, making the scheduling of IHLs challenging to handle. The optimal scheduling of BTs is confronted with the aforementioned computational complexity and uncertainty issues. First, the optimal scheduling of the on/off status of BTs is a complex integer programming problem. Second, since BTs are generally placed in the open air of an industrial site [12], the temperature change rate of a BT is easily influenced by weather conditions, such as temperature of outside ambiance. This could prevent the day-ahead scheduling from being implemented in reality.

Currently, exact methods and heuristic algorithms serve as the two primary categories of strategies to address the NP-hard problems. Many scholars have conducted extensive research on exact methods aiming at enhancing the solving efficiency of the NP-hard problems, including cutting plane-based methods [

13], branch- and bound-based methods [14], and methods based on Lagrange relaxation [15]. In addition, integer programming solvers, such as CPLEX and Gurobi, can automatically select and apply suitable exact methods to solve NP-hard problems, which are widely used in industry and academia [16]. However, while exact methods guarantee the global optimum, they may become impractical for large or complex problems due to their high computational demands. To solve the NP-hard problems in a more efficient manner, heuristic algorithms are generally employed. However, a major limitation of heuristic algorithms is that the optimality cannot be guaranteed [17]. Moreover, for all the aforementioned methods, solving efficiency decreases with the increasing scale of the NP-hard problems.

To address the calculation complexity arising from a large number of BTs and attain the optimal solution, a two-layer control model is established in this paper, which decomposes a control problem into two layers of sub-problems [

18]. The upper-layer control reduces the scale of the model by clustering the control objects to decrease the number of control variables. The lower-layer control distributes the upper-layer control commands to each control object. For instance, a two-layer clustering algorithm is developed in [19] to improve the computing efficiency of coordinated scheduling for electric vehicles in multi-microgrid systems. Reference [20] controls the electric devices by aggregating them with similar physical properties to provide ancillary services and promote supply-demand balance. Reference [21] introduces a clustering control method for controllable household load management, aiming at reducing peak consumption and achieving appropriate economic benefits. However, we have found no research focusing on the clustered analysis of IHLs, e.g., BTs, for expedited scheduling purposes. This is because the temperature transfer process of BTs is nonlinear over time, and the temperature distribution in the population of BTs is discrete, making the analysis of clustered temperature transfer characteristics challenging, especially in terms of average temperature constraint.

Moreover, there have been no studies considering the impact of external uncertainty factors on the temperature change rate of BTs in the scheduling processes. To cope with uncertainties, robust optimization is widely used, which is a mathematical method aiming to develop models that generate solutions that are less sensitive to the uncertainties [

22]. Previous studies typically describe uncertainty distributions using geometric shapes, such as box [23], ellipse [24], and polygon [25], which exhibit symmetrical spreading from the center. Therefore, they often encounter issues related to over-conservatism and reliance on a priori domain-specific knowledge. In real-world problems, uncertainty distributions can be extremely complex and challenging to be accurately captured through the traditional shapes. Therefore, applying a technique that can adaptively construct the uncertainty set while adjusting the level of conservatism is highly valuable. In this paper, support vector clustering (SVC) [26] is adopted to construct the uncertainty set due to its ability to capture uncertainties more flexibly based on the real historical data. Notably, the partial derivatives of the bitumen temperature with respect to the external uncertainty factors are derived to ensure that each element of the historical data samples has the same unit for the SVC implementation. Meanwhile, to address the computational burden caused by the soft margins of uncertain sets, a piecewise linear kernel (PLK) [27] is introduced in the SVC technique to approximate the uncertainty sets as convex polyhedrons. The implementation of the PLK-based SVC technique enables the flexible adjustment of the conservatism of the uncertainty set by modifying the regularization parameter. This adjustment is crucial for achieving the balance between the robustness and performance of the optimization results.

The main contributions of this paper are provided as follows.

1) The clustered temperature transfer characteristics of BTs are investigated with the aim of accelerating scheduling.

2) A two-layer deterministic scheduling (TDS) model is proposed to speed up the solution efficiency for the scheduling of BTs. In this model, the number of control variables is significantly reduced and remains unaffected by the increasing scale of BTs.

3) A data-driven PLK-based SVC technique is developed to deal with the uncertainties in the two-layer data-driven robust scheduling (TDRS) model of BTs. The shape of the uncertainty set can be adaptively constructed, and the level of its conservatism can be flexibly adjusted.

The remainder of the paper is arranged as follows. Section II introduces the modeling and control of BTs. Section III devises the TDS model for scheduling BTs. Section IV develops the two-layer model into a robust optimization format to deal with the uncertainties. Section V discusses the simulation results. Finally, Section VI concludes this paper and presents the future directions.

II. Modelling and Control of BTs

A. Industrial Sites

As shown in Fig. 1, an industrial site with BT, base load, and onsite photovoltaic (PV) generation can be treated as a microgrid connected to the power grid. In this industrial site, the participation of BTs in demand response is of great significance to accommodate onsite PV generation and reduce the impact of peak load on the power grid.

Fig. 1  An industrial site with BT, base load, and onsite PV generation.

B. Temperature Control of BTs

To ensure the availability of bitumen, it is necessary to keep it in BTs in a liquid state and, simultaneously, ensure that its temperature is within an allowable range [

12]. When the electrical switch of an industrial BT is turned on, the heat system is activated, causing the bitumen temperature to rise. Conversely, it naturally decreases when the electrical switch is turned off. Mathematically, the relationship between the heat transfer process of a BT and the on/off state (x) of its heater is given as:

Pnet=Pabsorb-Ploss=Pratex-UAT-Tamb (1)

In (1), x equals 1 if the heater is switched on and 0 if off. Meanwhile, Pnet decides the temperature change rate within the BT:

Pnet=cvmdTdt (2)

By combining (1) and (2), the temperature change rate of the BT can be represented by its switching state:

dTdt=Pratexcvm-UAT-Tambcvm (3)

Equation (3) shows that the temperature change rate is also related to the current temperature of bitumen, so it is constantly changing, even with the same switching state.

Figure 2 illustrates the temperature control process of a BT. The bitumen temperature increases when the heater is activated, while it decreases when the heater remains off.

Fig. 2  Temperature control process of a BT.

In Fig. 2, the electrical heater is set to be automatically turned off when the bitumen temperature touches the upper limit Tup at time t1, t 3, or t 5, while automatically turned on when the bitumen temperature reaches the lower limit Tdown at time t2. Moreover, the heater can also be manually switched on/off at any chosen time, e.g., t4.

III. Formulation of TDS Model

A. Original Deterministic Scheduling (ODS) Model

The ODS model is the conventional model for the scheduling of BTs. In this paper, the objective of the ODS model is to minimize the peak-to-valley difference of the electricity exchange between the industrial site and the external power grid.

minf=max1hHn=1Nxn,hPrate+PBLh-PPVh-min1hHn=1Nxn,hPrate+PBLh-PPVh (4)

The temperature transfer constraint of each BT and the temperature limit in each BT are given as:

Tn,h=Tn,h-1+Pratexn,h-UATn,h-1-TambcvmΔt    n,h (5)
TdownTn,hTup    n,h (6)

The ODS model presented above is an integer programming problem, with the number of control variables being N×H. Therefore, the solution time grows exponentially as the number of BTs increases.

B. TDS Model

To deal with the computational burden as the number of BTs increases, the TDS model is proposed.

1) Upper-layer Clustered Optimization

The upper-layer clustered optimization in the TDS model aims to decide the total number of BTs turned on at each time slot, by contrast to deciding the on/off status of each individual BT in the ODS model. Therefore, the number of decision variables is much reduced and does not change with the increasing numbers of BTs.

In this layer, the objective is formulated as:

minf=max1hHxhPrate+PBLh-PPVh-min1hHxhPrate+PBLh-PPVh (7)

Accordingly, the constraints of all the individual BTs are also replaced by those specifying the average temperature change of the whole BT population, i.e., the clustered temperature transfer process. Specifically, by summing (5) of each BT and calculating the mean, the constraint for the average temperature of the entire BT population is given as:

T¯h=T¯h-1+PratexhN-UAT¯h-1-TambcvmΔt    h (8)

The temperatures of BTs are actually discretely distributed around the average temperature. To ensure that the temperature of each BT does not go beyond the upper or lower limit, a certain gap needs to be kept between the average temperature of BTs and Tup and Tdown:

Tdown+ΔTT¯hTup-ΔT    h (9)

Remark: the control variables of this upper-layer clustered optimization model are xhh=1H. The total number of these control variables is H, which will not be influenced by the number of BTs. Therefore, compared with the ODS model, the calculation time can be significantly reduced when a large number of BTs are dealt with.

2) Lower-layer State Distribution

Based on the obtained upper-layer clustered optimization results, the lower-layer state distribution aims to decide which specific BTs are turned on at each time slot. Basically, it should be satisfied that the total number of BTs turned on after state distribution equals the obtained upper-layer clustered optimization results:

xh=n=1Nxn,h    h (10)

However, there are many feasible state distribution results under the satisfaction of (10). Therefore, this paper proposes a lower-layer state distribution principle, where BTs with lower temperatures are preferentially turned on at any time step. This ensures that the temperature of each BT remains as far from the temperature limits of Tup and Tdown as possible.

minf=1Nn=1NTn,h-T¯h2    hs.t.  (5), (6), (10) (11)

Mathematically, the state distribution principle is equivalent to minimizing the variance of temperatures of BTs during different time periods, as presented above.

C. Analysis on Gap Magnitude

In this subsection, a detailed analysis is provided to acquire the magnitude of the gap ΔT in (9). For a population of BTs, two arbitrary BTs are taken out for analysis. Based on the lower-layer state distribution, there are three combinations of the on/off states for the two BTs at any time step. ① Scenario 1: both BTs in the on state. ② Scenario 2: both BTs in the off state. ③ Scenario 3: the BT with a lower temperature in the on state while the BT with a higher temperature in the off state.

In Scenario 1, the BT with a lower temperature will heat up faster according to (5), as shown in Fig. 3(a). The red line represents the temperature transfer process of the BT with a higher initial temperature, while the blue line represents the temperature transfer process of the BT with a lower initial temperature. In Scenario 2, the BT with a higher temperature will cool down faster according to (5), as shown in Fig. 3(b). Hence, in both Scenarios 1 and 2, the temperature difference between the two BTs will narrow down after the subsequent time slot. In Scenario 3, assume that the absolute value of the initial temperature difference between the two BTs is ΔTinitial. If ΔTinitialΔT, the temperature difference will narrow down after the subsequent time slot, as shown in Fig. 3(c). If ΔTinitial<ΔT, the temperature difference may increase after the subsequent time slot, as shown in Fig. 3(d). From Fig. 3(d), it can be deduced that the maximum temperature difference between the two BTs after the subsequent time slot occurs when ΔTinitial0, as depicted in Fig. 3(e). The value of this maximum temperature difference between the two BTs is represented as ΔT. ΔTup is the temperature increase in the next Δt period for the BT with initial temperature T1,h, and ΔTdown is the temperature decrease for the BT with initial temperature T2,h.

Fig. 3  Temperature transfer process of two BTs. (a) Scenario 1. (b) Scenario 2. (c) Scenario 3 when ΔTinitialΔT. (d) Scenario 3 when ΔTinitial<ΔT. (e) Scenario 3 when ΔTinitial0.

ΔT=ΔTup+ΔTdown=Prate-UAT1,h-TambcvmΔt+         UAT2,h-TambcvmΔtPrateΔtcvm (12)

For the two BTs, T1,h is less than and infinitely approaches T2,h, i.e., T1,hT2,h. Hence, UAT2,h-T1,h/cvm0, and therefore ΔTPrateΔt/cvm.

Summarizing Scenarios 1-3, it can be deduced that if ΔTinitial<ΔT, the temperature difference between the two BTs after the subsequent time slot must be less than ΔT. Thus, as long as the very initial width of the temperature distribution of the BT population before the first time step is smaller than ΔT, it will not exceed ΔT at any time step. Therefore, the value of ΔT can be treated as the magnitude of the gap in (9). When constraint (9) is respected and the lower-layer state distribution is applied, the temperature of each individual BT will not go beyond the temperature limits throughout the scheduling horizon.

To conclude the formulation, the objective of the upper-layer clustered optimization is (7), which is subject to constraints (8), (9), and (12). After the upper-layer clustered optimization, the total number of on states at each time step can be calculated, and then the lower-layer state distribution is executed. The objective of lower-layer state distribution is (11), which is subject to constraints (5), (6), and (10).

IV. Formulation of TDRS Model

A. Uncertainty Analysis

The uncertainties of BTs lie in their clustered temperature transfer process. According to (8), there are two uncertain factors that influence the temperature change rate of BTs at an industrial site, specifically U and Tamb. For example, U would be higher in rainy and snowy weather while lower in sunny and dry weather. In addition, there are also forecasting errors for Tamb. Unlike the uncertainties of PV and base load, the uncertainties of BTs may cause automatic on/off switching of their electric heaters, leading to significant deviations in the actual clustered temperature transfer process of BTs from the day-ahead scheduling. This can result in fluctuations in the actual power consumption curve of BTs, deteriorating the execution results, which will be discussed in detail in Section V-B.

B. SVC-based Uncertainty Set Construction

To deal with the uncertainty of U and Tamb by using the robust optimal scheduling method, it is necessary to construct a convex two-dimensional uncertainty set for them.

1) Data-driven SVC Technique

In specific, assuming that M historical data samples umm=1M of the uncertain parameters considered are available, the objective of the SVC model is to seek a sphere that tries to enclose all the data samples with minimal volume and acts as the uncertainty set [

28]:

minR,ξR2+1Mvm=1Mξm (13)

s.t.

ϕum-P2R2+ξm    m (14)

In (13) and (14), ξmm=1M is adopted as the soft margins to accommodate outliers of historical data samples. v>0 is used to penalize outliers. The elements in um have the same unit. Moreover, the mapping function ϕ is designed to achieve the transformation from low-dimensional space to high-dimensional space, facilitating the enclosure of umm=1M.For solving the above SVC model, it is reformulated into a Lagrange function [

29] in (15), in which αm and βm are the mth Lagrange multipliers.

LP,R,ξ,α,β=R2+1Mvm=1Mξm-m=1MαmR2+ξm-ϕum-P2-m=1Mβmξm (15)

To solve (15), the following Karush-Kuhn-Tucker (KKT) conditions [

30] should be met:

LR=0m=1Mαm=1 (16)
LP=0P=m=1Mαmϕum (17)
Lξm=0αm+βm=1Mv (18)

Putting the KKT conditions into (15), the dual problem [

31] can be obtained as:

minαi=1Mj=1MαiαjϕuiTϕuj-i=1MαiϕuiTϕui (19)

s.t.

0αm1Mv    m (20)
m=1Mαm=1 (21)

Generally, to construct a convex two-dimensional uncertainty set, it is necessary to give the mapping function ϕ. However, in nonlinear classification problems, obtaining the appropriate mapping function from low-dimensional space to high-dimensional space may be challenging, and sometimes even infeasible [

32].

In this paper, the kernel function Κ,  is adopted to address this issue. Kernel function can solve the nonlinear classification problem in the low-dimensional space without explicitly giving the mapped function [

33]. Here, the PLK is selected as the kernel trick to facilitate solving the above dual problem, which implicitly calculates the inner product between data points in (19) as [27]:

Kui,uj=ϕuiTϕui=k=1Kminuik,ujk (22)

After calculating the dual problem, the related correlations of the results are described as [

26]:

ϕum-P2>R2αm=1Mv,βm=0ϕum-P2=R20<αm<1Mv,0<βm<1Mvϕum-P2<R2αm=0,βm=1Mvm (23)

where ϕum-P2>R2 indicates that um is an outlier outside the convex uncertainty set; ϕum-P2=R2 indicates that um is a boundary support vector, which forms the boundaries of the convex uncertainty set; and ϕum-P2<R2 indicates that um is a vector within the convex uncertainty set.

Therefore, the convex uncertainty set can be obtained as:

umϕum-P2R2,m (24)

In robust scheduling, due to the convexity of the constructed uncertainty set, it is only needed to ensure that the boundary support vectors satisfy the constraints.

2) Uncertainty Set Construction for Industrial Site

The above analysis is at the mathematical level. In actual industrial sites, the application of PLK-based SVC technique may encounter new issues. For example, as claimed in Section IV-B-1), um is a two-dimensional vector and its elements should have the same unit. However, this is not the case for U and Tamb in the scheduling of BTs.

To address this issue, this paper defines um consisting of the variation of the bitumen temperature due to the variation of U and Tamb, i.e., ΔTUm and ΔTTambm, which is calculated as:

ΔTUm=Um-UAT-TambΔtcvm (25)
ΔTTambm=Tambm-TambUAΔtcvm (26)

In (25) and (26), AT-TambΔtcvm and UAΔtcvm are the partial derivatives of the bitumen temperature with respect to U and Tamb, respectively, calculated from (5). Moreover, ΔTUm and ΔTTambm have the same unit (℃). Subsequently, the data-driven PLK-based SVC technique is used to generate the convex uncertainty set for scheduling BTs. Assuming the obtained boundary support vectors are donated as ΔTUq,ΔTTambqq=1Q, the corresponding historical data Uq,Tambqq=1Q can be used for the TDRS model.

C. TDRS Model

The objective function of the upper-layer clustered optimization is the same as (7). The consideration of uncertain factors is mainly reflected in the constraints, which are provided as:

T¯hq=T¯h-1q+PratexhN-UqAT¯h-1q-TambqcvmΔt    q,h (27)
Tdown+ΔTT¯hqTup-ΔT    q,h (28)

It should be noted that, the actual U and Tamb vary at different time, while the fixed parameter values are utilized in both the TDS model and each scenario of the TDRS model. If the temporal variation characteristics of Uq and Tambq are considered within each time slot, the dimension of the constructed uncertainty set should be 2×H, resulting in too many scenarios to be considered in the robust model and making the computation more complex. In fact, the fixed values of Uq and Tambq in the uncertainty set can also achieve satisfactory performance in addressing the uncertainties of BTs in the robust model, as will be demonstrated in Section V-D.

To conclude, firstly, the SVC-based uncertainty set construction is executed considering robustness. Secondly, based on the obtained boundary support vectors of the constructed uncertainty set, the TDRS model is implemented. The objective of the upper-layer clustered optimization is (7), which is subject to the constraints (12), (27), and (28). After obtaining the number of BTs turned on at each time slot considering uncertainties, the lower-layer state distribution is then executed. The objective of lower-layer state distribution is (11), subject to constraints (5), (6), and (10).

V. Case Study

In this section, the TDS model and the TDRS model are applied to the industrial system introduced in Section II-A. All numerical simulations are conducted on a laptop equipped with a 2.60 GHz i7 CPU processor and 8 GB RAM. The CPLEX solver is executed in MATLAB to solve the models. The parameter values used in the cases are listed in Table I, sourced from real data provided by the KVM UK Ltd [

4].

TABLE I  Parameter Values
ParameterValueParameterValue
U (kW·m-2·K-1) 7.75×10-3 Prate (kW) 120
cv (kJ·kg-1·K-1) 1.34 Tdown (℃) 150
Tup (℃) 180 m (kg) 21500
M 400 Δt (s) 900
A (m2) 36 H 96

A. Results of TDS Model

In this subsection, different existing methods for solving the ODS model, such as the CPLEX solver (with the optimality gap set at 0.01) and the PSO algorithm (with acceleration coefficients c1=1.8, c2=2.2, and inertia weight w=0.3), are compared with the proposed TDS model, which is also solved by the CPLEX solver with the same optimality gap. Moreover, the penalty function method is utilized in the PSO algorithm to handle the constraints.

Table II presents the comparison results of different optimal scheduling methods for N=10, 15, and 20, respectively. In terms of calculation time, the CPLEX solver outperforms the PSO algorithm when solving the ODS model. Nevertheless, when N=20, the calculation time of the CPLEX solver exceeds 2 hours, which is difficult to execute in actual industrial processes, e.g., in a receding horizon control manner. The calculation time of the TDS model indicates that it not only has a much faster computation speed but also remains unaffected by the increasing scale of BTs.

TABLE II  Comparison Results of Different Optimal Scheduling Methods
MethodCalculation time (s)Peak-to-valley difference (MW)
N=10N=15N=20N=10N=15N=20
CPLEX solver 214.38 968.47 7359.65 0.3872 0.1552 0.1136
PSO algorithm 6627.23 18223.71 0.4335 0.1871
TDS model 13.27 13.92 13.46 0.4221 0.1934 0.1136

The TDS results and the corresponding temperature transfer process of BTs when N=20 are provided, as shown in Fig. 4. Figure 4(a) depicts the total consumption curve of 20 BTs, which is approximately consistent with the curve of net PV generation (PV generation subtracted by base load). Additionally, the curve of the power exchange with the main grid is almost flat, achieving the scheduling objective of minimizing the peak-to-valley difference. Figure 4(b) exhibits the temperature transfer process of BTs. Each colored curve in Fig. 4(b) represents the temperature fluctuation of a BT. Meanwhile, the average temperature curve (represented by the black line) and the magnitude of the gap ΔT are also exhibited in Fig. 4(b). During approximately 10:00 and 16:30, the average temperature of BTs is increasing, which coincides with a relatively higher power consumption of BTs during this period. On the contrary, during other periods, the decreasing average temperature of BTs coincides with relatively lower power consumption values. ΔD0 in Fig. 4(b) shows the difference of the maximum and minimum temperatures of BTs during the entire temperature transfer process. Table II shows that, when N=20, the TDS achieves the same peak-to-valley difference as the ODS model. However, when N=10 and N=15, the TDS model achieves inferior peak-to-valley differences compared with the ODS model. Moreover, as indicated in Table II, the peak-to-valley difference increases as N decreases. The reasons are analyzed below.

Fig. 4  TDS results and corresponding temperature transfer process of BTs when N=20. (a) Total consumption of 20 BTs. (b) Corresponding temperature transfer process.

Figure 5 shows the TDS results and the corresponding temperature transfer process of BTs when N=10. Due to the limited scheduling capacity of 10 BTs, the curve of the power exchange with the main grid illustrated in Fig. 5(a)) exhibits a larger peak-to-valley difference than that in Fig. 4(a). Moreover, according to (9), the average temperature curve of BTs in the TDS results is confined between the dashed red lines. This leads to the temperature variation curve of each BT not being able to fall within the region of dashed black lines in Fig. 5(b). Thus, the temperature variation space of BTs in the TDS model is a little smaller than that in the ODS model, resulting in an inferior result in the peak-to-valley difference. Nevertheless, the gap in peak-to-valley difference between these two models is only 0.0349 MW, which is certainly industrially acceptable given the significant drop in computing duration. Similar analysis can be applied to the TDS results when N=15, which is not provided here.

Fig. 5  TDS results and corresponding temperature transfer process of BTs when N=10. (a) Total consumption of 10 BTs. (b) Corresponding temperature transfer process.

B. Influence of Uncertain Factors

Due to the influence of uncertain factors, the control commands of BTs obtained from the TDS model may lead to constraint violation. Figure 6 shows the empirical U and the forecasted Tamb on a rainy day, along with their respective actual values. On this day, the actual values of U are higher than the empirical value due to the wet weather.

Fig. 6  Values of uncertain factors on a rainy day. (a) U. (b) Tamb.

Figure 7(a) illustrates the actual total consumption of 20 BTs on this day, which are based on the control commands obtained from the TDS model. Figure 7(b) displays the corresponding temperature transfer process. In the actual scenario shown in Fig. 6, due to the higher U and lower Tamb before approximately 09:00, the heat of BTs dissipates faster than expected. Therefore, the actual temperatures of BTs drop more rapidly compared with the TDS results illustrated in Fig. 4(b), leading to the temperatures of BTs approaching Tdown at about 09:00, as shown in Fig. 7(b). The electrical heater of some BTs automatically turns on when their temperatures reach Tdown, leading to obvious fluctuations during the following period in the actual power consumption curve of BTs, as circled in Fig. 7(a).

Fig. 7  Actual execution results of 20 BTs obtained by TDS model under uncertainty. (a) Actual total consumption of 20 BTs. (b) Corresponding temperature transfer process.

C. Generation of Uncertainty Sets

In this paper, 400 historical data samples of um are collected from the KVM UK Ltd [

4] and shown in Fig. 8. From Fig. 8, it can be observed that the uncertainty in U may cause a maximum deviation in the bitumen temperature of about 0.188 ℃ per Δt=15 min, while the uncertainty in Tamb may result in a maximum deviation of about 0.034 ℃ per Δt. This implies that compared with the uncertainty in Tamb, the uncertainty in U is likely to cause more deviations in the bitumen temperature. If these uncertain factors are not effectively considered in the scheduling models, they could lead to accumulated temperature deviations in actual execution, as illustrated in Fig. 7(b).

Fig. 8  Uncertainty sets of um. (a) v=0.25. (b) v=0.125.

For example, on a rainy day, the rain may last for several hours, e.g., 6 hours. During this rainy period, the higher U can lead to the accumulated temperature deviation reaching (0.188×24) =4.512 ℃. Considering that the allowable temperature range of BTs is from 150 ℃ to 180 ℃, this accumulated temperature deviation is notable and may result in the bitumen temperature exceeding the allowable range. If that happens, the BTs, detecting this deviation via temperature sensors, would have to move away from ODS to correct the bitumen temperature back to the allowable range, thus leading to larger power fluctuations of BTs, as analyzed in Section V-B.

Figure 8 illustrates the generated uncertainty sets by using the PLK-based SVC technique. The conservatism of the uncertainty set can be increased by reducing the regularization parameter v, as introduced in (13). Through utilizing the boundary support vectors in the TDRS model, the influence of uncertain factors in the actual execution of control commands of BTs can be considered.

D. Results of TDRS Model

The actually executed power consumption curve of 20 BTs based on the TDRS model with v=0.25 is illustrated in Fig. 9(a). By considering the uncertainty set in the TDRS model, no BT is unexpectedly turned on or off during the actual temperature transfer process, as shown in Fig. 9(b).

Fig. 9  Actual execution results of 20 BTs based on TDRS model with v=0.25. (a) Actual total consumption of 20 BTs. (b) Corresponding temperature transfer process.

It is noteworthy that in the TDRS model, due to the consideration of the uncertainty set, the temperature difference, ΔD3 (shown in Fig. 9(b)), is reduced compared with ΔD0(shown in Fig. 4(b)). This is because, to avoid violating the temperature limits in actual scenarios, the TDRS results turn on more BTs before approximately 10:00, and more BTs are turned off between approximately 10:00 and 16:30.

Figure 10(a) presents the actual total consumption results of 20 BTs based on the TDRS model with v=0.125. In this case, the conservatism of the uncertainty set is further increased compared with that with v=0.25, leading to a further reduction in the temperature difference ΔD4, compared with ΔD3.

Fig. 10  Actual execution results of 20 BTs based on TDRS model with v=0.125. (a) Actual total consumption of 20 BTs. (b) Corresponding temperature transfer process.

Consequently, more BTs are turned on before approximately 10:00, and more are turned off between approximately 10:00 and 16:30. As a result, the actual electricity exchange with the power grid obtained by the TDRS model with v=0.125 is higher than that with v=0.25 before approximately 10:00 and is lower than that between approximately 10:00 and 16:30.

Table III compares the actual peak-to-valley difference of TDS model and TDRS model under different levels of conservatism. When the ellipse is treated as the uncertainty set, the least squares method is employed to form an elliptical region that covers all historical data samples. Moreover, when the box is treated as the uncertainty set, the maximum and minimum values of umm=1M in each dimension are considered as the boundaries of the box. Due to the traditional uncertainty sets (ellipse and box) covering all historical data samples, they are more conservative than the uncertainty sets generated by using the PLK-based SVC technique. Therefore, initially, with an increase in conservatism (decrease in v), the actual peak-to-valley difference decreases. However, as conservatism further increases, e.g., when v=0.125 or when the ellipse or box is used as the uncertainty set, the peak-to-valley characteristic gradually deteriorates due to over-conservatism. The comparison results show that when v=0.3, it is an appropriate level of conservatism that enables the proposed TDRS model to achieve a better actual peak-to-valley difference.

TABLE III  Comparison of Actual Peak-to-valley Difference
ModelPeak-to-valley difference (MW)
TDS model 0.3820
TDRS model with v=0.4 0.2181
TDRS model with v=0.35 0.1847
TDRS model with v=0.3 0.1695
TDRS model with v=0.25 0.1741
TDRS model with v=0.2 0.2208
TDRS model with v=0.125 0.2914
TDRS model (with an ellipse as uncertainty set) 0.4003
TDRS model (with a box as uncertainty set) 0.4198

Remark: in Table III, the TDRS model with different values of v is executed in a specific actual scenario depicted in Fig. 6. To obtain the most suitable level of conservatism for the TDRS model, it is needed to examine its performance across various actual scenarios.

In this paper, 100 historical data samples are randomly selected from the 400 samples depicted in Fig. 8. Then, the performance of the TDRS model with different values of v is examined for various actual scenarios corresponding to the selected samples. In general, a robust model with lower conservatism exhibits superior performance in the scenarios with lower uncertainty; conversely, it performs better in the scenarios with higher uncertainty when featuring higher conservatism. Figure 11 illustrates that the TDRS model with v=0.3 outperforms the model with other values of v in 35 out of 100 actual scenarios, making it the recommended choice in this paper.

Fig. 11  Number of actual scenarios achieving a better peak-to-valley difference under TDRS model with different values of v.

VI. Conclusion

In this paper, a two-layer TDRS model is proposed for tapping the flexibility contained in IHLs. To deal with the calculation difficulty and uncertainties, the clustered temperature transfer process of BTs is investigated, and a PLK-based SVC technique is utilized to deal with the uncertainties in robust optimal scheduling. Compared with the ODS model, the proposed TDRS model can significantly reduce the computation time of scheduling 20 BTs from over 2 hours to 13.46 s, with both models being solved using the CPLEX solver. Moreover, by comparing the results of the TDRS model in various actual scenarios, a recommended level of conservatism (v=0.3) is adopted to mitigate the influence of uncertain factors in the actual execution of day-ahead scheduling.

The appropriate conservatism of the TDRS model varies across various actual scenarios, as shown in Fig. 11, and the characteristics of actual scenarios are correlated with external weather conditions. Therefore, combining the local weather monitoring system to establish an adaptive conservatism adjustment strategy may further enhance the effectiveness of the TDRS model, which is considered a direction for future research. Moreover, given that BTs are equipped with real-time temperature sensors, it is possible to iteratively correct the accumulated temperature deviation of BTs within a receding-horizon framework, e.g., real-time predictive control. This presents another promising direction for future research.

Nomenclature

Symbol —— Definition
A. —— Parameters
Δt —— Length of each time slot
ΔT —— Magnitude of gap
ΔTup —— Temperature increase during a Δt period for a bitumen tank (BT)
ΔTdown —— Temperature decrease during a Δt period for a BT
ΔTUm —— Temperature transfer deviation of BTs caused by uncertainty of Um
ΔTTambm —— Temperature transfer deviation of BTs caused by uncertainty of Tambm
α, β —— Lagrange multipliers
ξ —— Slack variable
ξm —— The mth slack variable
A —— Area of a BT
cv —— Heat capacity of bitumen
H —— Number of time slots in time horizon
K —— Dimension of historical data sample
M —— Number of historical data samples
m —— Mass of bitumen in a BT
N —— Number of BTs
Pabsorb —— Heat absorption rate
Ploss —— Heat loss rate
Pnet —— Heat transfer rate
Prate —— Rated heat power
P —— Center of sphere
Q —— Number of boundary support vectors
R —— Radius of sphere
T —— Bitumen temperature
Tamb —— Temperature of outside ambiance
Tup —— Upper limit of bitumen temperature
Tdown —— Lower limit of bitumen temperature
Tambm —— The mth historical data of Tamb
U —— Overall heat transfer coefficient
Um —— The mth historical data of U
um, ui, uj —— The mth, ith, jth historical parameter samples
uik —— The kth element of ui
ujk —— The kth element of uj
v —— Regularization parameter
B. —— State Variables
PBLh —— Base load of industrial site at time h
PPVh —— Power generation of local photovoltaic (PV) at time h
T¯h —— Average temperature of BTs at time h
Tn,h —— Temperature of the nth BT at time h
T¯hq —— Average temperature of BTs corresponding to the qth boundary support vector at time h
C. —— Decision Variables
x —— On/off state of a BT
xn,h —— State of the nth BT at time h
xh —— Total number of BTs turned on at time h

References

1

R. Hanna and D. G. Victor, “Marking the decarbonization revolutions,” Nature Energy, vol. 6, no. 6, pp. 568-571, Jun. 2021. [Baidu Scholar] 

2

C. Dang, J. Zhang, C.-P. Kwong et al., “Demand side load management for big industrial energy users under blockchain-based peer-to-peer electricity market,” IEEE Transactions on Smart Grid, vol. 10, no. 6, pp. 6426-6435, Nov. 2019. [Baidu Scholar] 

3

Precedence Research. (2023, Jan.). Bitumen market size, 2022 to 2032. [Online]. Available: https://www.precedenceresearch.com/bitumen-market [Baidu Scholar] 

4

Asphalt Industry Alliance. (2023, Dec.). Market changes reflected in latest UK asphalt industry report. [Online]. Available: https://www.asphaltuk.org/key-facts/ [Baidu Scholar] 

5

TEC Container Solutions. (2022, Dec.). How is bitumen heated? [Online]. Available: https://teccontainersolutions.com/2022/12/how-is-bitumen-heated/ [Baidu Scholar] 

6

MBA Plant. (2022, Nov.). MBA asphalt plants. [Online]. Available: https://mbaplant.com/en/asphalt-plant/ [Baidu Scholar] 

7

X. Xiao, D. Cai, L. Lou et al., “Application of asphalt based materials in railway systems: a review,” Construction and Building Materials, vol. 304, p. 124630, Oct. 2021. [Baidu Scholar] 

8

Y. Zhou, M. Cheng, and J. Wu, “Enhanced frequency response from industrial heating loads for electric power systems,” IEEE Transactions on Industrial Informatics, vol. 15, no. 6, pp. 3388-3399, Jun. 2018. [Baidu Scholar] 

9

M. Cheng, J. Wu, S. J. Galsworthy et al., “Power system frequency response from the control of bitumen tanks,” IEEE Transactions on Power Systems, vol. 31, no. 3, pp. 1769-1778, May 2015. [Baidu Scholar] 

10

C. Wu, H. Han, S. Gao et al., “Coordinated scheduling for multi-microgrid systems considering mobile energy storage characteristics of electric vehicles,” IEEE Transactions on Transportation Electrification, vol. 9, no. 1, pp. 1775-1783. Mar. 2022. [Baidu Scholar] 

11

F. Verastegui, A. Lorca, D. E. Olivares et al., “An adaptive robust optimization model for power systems planning with operational uncertainty,” IEEE Transactions on Power Systems, vol. 34, no. 6, pp. 4606-4616, Nov. 2019. [Baidu Scholar] 

12

A. Hamed, E. Shaban, R. Darwish et al., “Design and implementation of discrete PID control applied to Bitumen tank based on new approach of pole placement technique,” International Journal of Dynamics and Control, vol. 5, pp. 604-613, Aug. 2017. [Baidu Scholar] 

13

A. Testa, A. Rucco, and G. Notarstefano, “Distributed mixed-integer linear programming via cut generation and constraint exchange,” IEEE Transactions on Automatic Control, vol. 65, no. 4, pp. 1456-1467, Apr. 2019. [Baidu Scholar] 

14

F. Najafi and M. Fripp, “Stochastic optimization of comfort-centered model of electrical water heater using mixed integer linear programming,” Sustainable Energy Technologies and Assessments, vol. 42, p. 100834, Dec. 2020. [Baidu Scholar] 

15

M. A. Bragin and E. L. Tucker, “Surrogate “level-based” Lagrangian relaxation for mixed-integer linear programming,” Scientific Reports, vol. 12, no. 1, p. 22417, Dec. 2022. [Baidu Scholar] 

16

C. Wu, S. Jiang, S. Gao et al., “Event-triggered model predictive control for dynamic energy management of electric vehicles in microgrids,” Journal of Cleaner Production, vol. 368, p. 133175, Sept. 2022. [Baidu Scholar] 

17

K. Amine, “Multiobjective simulated annealing: principles and algorithm variants,” Advances in Operations Research, vol. 2019, p. 8134674, May 2019. [Baidu Scholar] 

18

M. Elkazaz, M. Sumner, E. Naghiyev et al., “A hierarchical two-stage energy management for a home microgrid using model predictive and real-time controllers,” Applied Energy, vol. 269, p. 115118, Jul. 2020. [Baidu Scholar] 

19

C. Wu, H. Han, S. Gao et al., “Coordinated scheduling for multimicrogrid systems considering mobile energy storage characteristics of electric vehicles,” IEEE Transactions on Transportation Electrification, vol. 9, no. 1, pp. 1775-1783, Mar. 2022. [Baidu Scholar] 

20

X. Wu, L. You, R. Wu et al., “Management and control of load clusters for ancillary services using internet of electric loads based on cloud-edge-end distributed computing,” IEEE Internet of Things Journal, vol. 9, no. 19, pp. 18267-18279, Oct. 2022. [Baidu Scholar] 

21

R. Pereira, A. Fagundes, R. Melício et al., “A fuzzy clustering approach to a demand response model,” International Journal of Electrical Power & Energy Systems, vol. 81, pp. 184-192, Oct. 2016. [Baidu Scholar] 

22

B. Zhao, T. Qian, W. Tang et al., “A data-enhanced distributionally robust optimization method for economic dispatch of integrated electricity and natural gas systems with wind uncertainty,” Energy, vol. 243, p. 123113, Mar. 2022. [Baidu Scholar] 

23

M. Zhang, J. Fang, X. Ai et al., “Partition-combine uncertainty set for robust unit commitment,” IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 3266-3269, Jul. 2020. [Baidu Scholar] 

24

H. Daneshvari and R. Shafaei, “A new correlated polyhedral uncertainty set for robust optimization,” Computers & Industrial Engineering, vol. 157, p. 107346, Jul. 2021. [Baidu Scholar] 

25

D. Bertsimas and M. Sim, “The price of robustness,” Operations Research, vol. 52, no. 1, pp. 35-53, Feb. 2004. [Baidu Scholar] 

26

A. Ben-Hur, D. Horn, H. T. Siegelmann et al., “Support vector clustering,” Journal of Machine Learning Research, vol. 2, pp. 125-137, Dec. 2001. [Baidu Scholar] 

27

P. Qiu, “A jump-preserving curve fitting procedure based on local piecewise-linear kernel estimation,” Journal of Nonparametric Statistics, vol. 15, no. 4-5, pp. 437-453, Oct. 2003. [Baidu Scholar] 

28

C. Shang, X. Huang, and F. You, “Data-driven robust optimization based on kernel learning,” Computers & Chemical Engineering, vol. 106, pp. 464-479, Nov. 2017. [Baidu Scholar] 

29

S. Liang, X. Zeng, and Y. Hong, “Distributed nonsmooth optimization with coupled inequality constraints via modified Lagrangian function,” IEEE Transactions on Automatic Control, vol. 63, no. 6, pp. 1753-1759, Jun. 2017. [Baidu Scholar] 

30

A. Sinha, T. Soun, and K. Deb, “Using Karush-Kuhn-Tucker proximity measure for solving bilevel optimization problems,” Swarm and evolutionary computation, vol. 44, pp. 496-510, Feb. 2019. [Baidu Scholar] 

31

H. Zhu, F. Ye, and E. Zhou, “Solving the dual problems of dynamic programs via regression,” IEEE Transactions on Automatic Control, vol. 63, no. 5, pp. 1340-1355, May 2017. [Baidu Scholar] 

32

D. A. Pisner and D. M. Schnyer, “Support vector machine,” in Machine Learning. Amsterdam: Elsevier, 2020, pp. 101-121. [Baidu Scholar] 

33

M. Shao, X. Wang, Z. Bu et al., “Prediction of energy consumption in hotel buildings via support vector machines,” Sustainable Cities and Society, vol. 57, p. 102128, Jun. 2020. [Baidu Scholar]