Next Article in Journal
Enhancing Linux System Security: A Kernel-Based Approach to Fileless Malware Detection and Mitigation
Previous Article in Journal
Blockchain Forensics: A Systematic Literature Review of Techniques, Applications, Challenges, and Future Directions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bi-Level Optimal Configuration of Electric Thermal Storage Boilers in Thermal–Electrical Integrated Energy System

School of Automation and Electrical Engineering, Inner Mongolia University of Science and Technology, Baotou 014010, China
*
Author to whom correspondence should be addressed.
Electronics 2024, 13(17), 3567; https://doi.org/10.3390/electronics13173567
Submission received: 24 July 2024 / Revised: 2 September 2024 / Accepted: 5 September 2024 / Published: 8 September 2024

Abstract

:
Electric thermal storage boilers (ETSBs) are important devices in enhancing the electric–thermal decoupling ability and spatiotemporal transfer of integrated energy system (IES), which is beneficial for improving system flexibility and energy utilization efficiency. In order to obtain more accurate and comprehensive results, a bi-level optimal model is proposed to study the site selection and capacity configuration of ETSB in IES based on the established mathematical model of ETSB. The objective of upper-level optimization of the model is obtaining the lowest energy supply cost when configuring the location and capacity of ETSB, while the lower-level model optimizes the operation scheduling with the goal of obtaining the lowest operational cost. The mixed-integer linear programming method and the genetic algorithm method are selected to obtain the optimal model. To illustrate the effectiveness and advantages of the proposed method, case studies are carried out. The optimal configuration scheme for an ETSB is obtained by comparing the lowest energy supply cost under different configuration parameters. Furthermore, the impact of an ETSB on the system is also analyzed based on the variations in energy balance, abandoned energy, and energy allocation before and after configuring the ETSB.

1. Introduction

To achieve its carbon peak and neutrality targets, the active decarbonization of energy systems has become inevitable in China. To this end, it is vital to promote shifting the primary energy source from fossil energy to renewable energy. Integrated energy systems, which seek to realize integrated energy production, distribution, and supply, are considered as a critical technology for promoting the consumption of renewable energy and have received widespread attention in the international energy field. However, it is nevertheless clear that the insufficient flexibility of IES caused by the imbalance of source–load matching and the volatility of renewable energy limits the development of IES [1]. To address this issue, IES must be able to further enhance the system’s adaptability and responsiveness to ensure the reliability of the energy supply.
Utilizing multi-energy coupling to achieve the conversion and complementarity of different types of energy is one of the key technologies used in IES to improve system flexibility and energy efficiency. Many scholars have conducted studies on this technology. Li et al. [2] built a simulation platform based on TRNSYS and Genopt software to deal with the problem of the coupling relationship between various sections of a system during the optimization of the regional IES capacity configuration. Lu et al. [3] established a correlation model of combined cooling, heating, and power coupled multi-energy system for optimization based on a bi-level model construction method. In ref. [4], a two-stage mixed-integer linear programming approach for district level multi-energy system planning was proposed to solve the problem of distributed renewable energy integration. In ref. [5], the installation configuration of gas turbines was determined using a multi-objective optimization model based on an energy hub by comparing the results regarding the economic, environmental, and energy consumption benefits simultaneously. The electric boiler (EB), as a typical electric thermal energy coupling device, is one of the main pieces of equipment used to decouple electric heat and clean energy heating. In recent years, various scholars have conducted research in this field. Meibom et al. [6] analyzed EB’s operating characteristics and the economic benefits of introducing EB in three district heating systems within the North European power system. The results show that this has significant beneficial effects on the wind power consumption capacity and system economy. Zhao et al. [7] presented an EB configuration scheme for increasing the electrical and thermal coupling ability, which could further improve the flexibility of combined heat and power plant (CHP) with wind power generation systems.
In addition, energy storage technology is another effective way for IES to improve the stability of energy supply systems, efficiently achieving the spatiotemporal migration of energy by storing and releasing it. In order to enrich the understanding of the effects of energy storage technology, Wang et al. [8] summarized the characteristics of multiple energy storage technologies and evaluated their ability to alleviate fluctuations and uncertainties in renewable energy. Bazdar et al. [9] presented a comprehensive review of technological features in compressed air energy storage and analyzed the latter’s role in micro-grids from the perspectives of integration potential, optimization design, and scheduling. Datta et al. [10] summarized and discussed the various applications of battery energy storage system technology, from power conversion systems to hybrid system configurations, in reducing the adverse impacts of renewable energy. Thermal energy storage enables the realization of the sustainable and economic operation of electric–heat integrated energy systems (EH-IES). To study the model theory and devise a configuration method for thermal storage system (TSS), Ren et al. [11] proposed an integrated optimization framework. Wei et al. [12] presented a two-stage optimal configuration model for TSS in IES based on analyzing the complex operational features.
In real projects, ETSB is always used to enhance the flexibility of the system with the dual function of EB and energy storage technology, converting surplus electricity into thermal energy and storing it for use in district heating. This is considered to be very promising for improving the flexibility of electric–thermal energy systems and suppressing system energy volatility [13]. Rakesh et al. [14] proposed integrating ETSB into the Danish energy network to increase the flexibility of the demand side of the system and simulated ETSB model in low-pressure residential buildings to evaluate the utilization potential. Zhao et al. [15] conducted a techno-economic analysis by comparing the schemes of heat-only boilers and ETSB for enhancing the flexibility of CHP, and their results indicated that the net annual revenue of the latter scheme was significantly higher. Based on the active wind curtailment model of ETSB, Lei et al. [16] proposed a strategy of wind curtailment and production consumption for a combined heating system using ETSB and a coal-fired boiler. Liu et al. [17] introduced an optimal economic dispatching model to minimize the operational cost and improve the energy utilization efficiency of CHPs by using ETSB.
The studies mentioned above have laid the theoretical and model foundations for the research conducted in this article. However, there are still some issues left unresolved.
(1) It is widely recognized that reasonable planning results are the foundation for achieving efficient system operation. However, almost all scholars just focus on the operational effect of ETSB on IES and a limited number of studies have been conducted on scientific planning and optimal configuration.
(2) ETSB includes three parameters: the electrical input power, the thermal output power, and the thermal storage capacity. Due to the diversity and coupling of energy supply methods, the operation and analysis of energy flows are complex in IES, meaning that all three parameters can affect the distribution of the system energy flow in IES during operation. However, existing studies usually considered only one or two ETSB parameters to analyze their effect on IES, which will inevitably lead to deviations in the analyzed results.
In response to the above problems, a bi-level optimal ETSB configuration method for electric–thermal IES is proposed in this paper by establishing a refined mathematics model of ETSB. We analyze the impact of ETSB on the EH-IES, which can enrich the theoretical study of IES, providing a more reasonable configuration scheme.
The remainder of this paper is organized as follows: the EH-IES model is presented and the refined mathematics model of ETSB is described in Section 2. Section 3 elaborates on the bi-level optimal configuration of ETSB. In Section 4, a case study of the optimal configuration is introduced to validate the effectiveness of the model and analyze the characteristics of the optimization results.

2. Model of IES

2.1. Typical Physical Architecture of IES

Through the coupling between different energy sources, IES is an integrated system of production and supply of various forms of energy, proposed to improve the efficiency and flexibility and increase the consumption of renewable energy. EH-IES can provide electricity and thermal energy to end users through energy conversion and storage, as shown in Figure 1.
As can be seen, some typical energy devices are contained in this system. CHP is adopted as the main unit of power and heat supply, and wind power (WP) and photovoltaic (PV) are included as the renewable energy units. These energy devices generate electricity and thermal energy, which are injected into the power and heating network, respectively, to meet the demands of power and thermal loads. To suppress the volatility of the source and load of the system, ETSB is introduced to improve the electric–thermal conversion.

2.2. Mathematical Model of ETSB

The energy conversion process of ETSB is shown in Figure 2. ETSB uses resistance or electromagnetic heating technology to convert excess electrical energy from the power network into thermal energy. In addition to meeting the needs of the heating network, the surplus thermal energy is stored in the TSS. Excess thermal energy from the heating network can be directly stored in the TSS. Thus, the ETSB can be regarded as an organic combination of EB and TSS units, which can achieve electric–thermal decoupling and spatiotemporal transfer [18]. Based on the operation principle of ETSB, the energy conversion process follows the following formulas.
Q E , t + P E B , t η E = Q s t o r e , t + Q o u t , t
where Q E , t is the input surplus thermal energy of the ETSB at time t, P E B , t represents the input power of the ETSB at time t, η E is the conversion efficiency of electricity to thermal energy, Q o u t , t is the thermal energy discharged by the ETSB at time t, and Q s t o r e , t is the thermal energy stored by the ETSB at time t, which can be determined by P E , t and Q o u t , t .
Moreover, the thermal storage state of ETSB can be expressed as Equation (2). It is easy to obtain that Q o u t , t is constrained by the thermal storage capacity of ETSB Q s t , c a .
Q s t o r e , t = Q s t o r e , t 1 ( 1 η s ) Q o u t , t + P E B , t η E + Q E , t
where Q s t o r e , t 1 and Q s t o r e , t 1 are the thermal storage states of ETSB at time t and t − 1, while η s is the thermal charge efficiency of ETSB.

3. Optimal Configuration Model

It is widely recognized that IES will play a key role in future energy development, and the planning and configuration of IES is important to optimize system operation. In this paper, the ETSB sites are selected and their capacity is configured using an existing EH-IES. These processes are typically carried out during expansion planning.
According to the mathematical description of ETSB mentioned above, it can be seen that ETSB has dual attributes relating to electricity and thermal energy. Without adding system lines or pipelines, the ETSB needs to be configured at the electricity–thermal coupling node in the IES, that is, at CHP nodes in the IES framework shown in Figure 1.
In the power subsystem, the ETSB acts as a load, while in the heating subsystem, it is equivalent to a thermal source. Therefore, configuring the ETSB in the existing IES system can provide the ability to transfer electricity to thermal energy, resulting in the electric–thermal decoupling effect. When the ETSB supplies thermal energy to the heating subsystem, consuming some surplus electricity, the energy flows will redistribute in the system based on the optimal operation. However, the configuration of the ETSB will also introduce certain system investment costs. Based on this, a bi-level optimal configuration model is proposed to plan and configure ETSB in EH-IES. The optimal scheme is shown in Figure 3.
In the upper level, the location selection and capacity configuration of the ETSB can be determined on the basis of the equipment and energy constraints with the goal of minimizing the energy supply cost of the system. Then, in the lower level, the energy flow distribution results of the system can be solved, and the operating status of each node can be obtained with the goal of obtaining the lowest operational cost. The results of the upper level will generally affect the objective function and constraints of the lower level, while the operation results are returned to the upper level, so as to achieve the interaction between the two levels of decisions. Figure 3 presents the bi-level optimization logic diagram.

3.1. Capacity Configuration Model of the Upper Level

3.1.1. Objective Function

The configuration of the ETSB in the EH-IES is optimized from the perspective of achieving optimal economy. The minimum total daily energy supply cost of the system is taken as the objective function in this study:
min C t o t a l = C i n v + C o p
where C i n v is the daily cost of ETSB investment and C o p is the operational cost of the system.
The daily investment cost of the ETSB is obtained by averaging the total investment cost of the ETSB for each operational day by considering the time value of the ETSB.
C i n v = C E T S B C R F n / n y
where C E T S B is the initial investment cost of the ETSB; n y is the annual number of operational days of the equipment; and C R F n is the capital recovery factor, which can be calculated as follows [19]:
C R F n = r ( 1 + r ) n ( 1 + r ) n 1
where r is the discount rate, set as 5% in this paper, while n is the lifetime of the ETSB.
The operational cost of the system C o p mainly includes the energy purchase cost C o p , p r and the operational maintenance cost C o p , o m . It was assumed that the IES runs in island mode in this paper, so the energy purchase cost only includes gas purchases. Moreover, the operational maintenance price of devices is calculated as follows:
C o p = C o p , p r + C o p , o m
C o p , p r = z Ω CHP γ g P g , p r , z
C o p , o m = a Ω CHP ε CHP P CHP , a + b Ω WP ε WP P WP , b + c Ω P V ε PV P PV , c + d Ω EBST ε ETSB P ETSB , d
where P CHP , P WP , P PV , and P ETSB , respectively, represent the power of the CHP, WP, PV, and ETSB; Ω CHP , Ω WP , Ω PV and Ω ETSB , respectively, represent the device collection of the CHP, WP, PV, and ETSB; γ g represents the gas purchase price; P g , p r , z represents the gas purchase volume of CHP z; and ε CHP , ε WP , ε PV , and ε ETSB , respectively, represent the operational maintenance cost related to the power of the device. Generally, the operational maintenance cost of the devices is mainly based on routine maintenance and annual overhaul. In this paper, the operational maintenance cost prices of the devices were determined by establishing the relationship between annual maintenance cost and output based on empirical data.

3.1.2. Constraint Function

Based on the multiple uncertainties regarding wind power and photovoltaic output, the constraints of the optimal configuration model of the system are as follows:
(1) Power balance constraints
Power balance refers to the power supply meeting the power demand of the IES. Thus, the power balance constraints of the system can be described as follows:
P CHP , t + P WP , t + P PV , t P EB , t = P L , t + P P , l o s s , t
where P CHP , t is the electrical power generated by the CHP at time t; P EB , t is the electrical power of the ETSB at time t; P WP , t and P PV , t represent the power output of the WP and PV at time t; P L , t represents the electrical load of the system at time t; and P P , l o s s , t represents the electrical loss of the system at time t.
(2) Thermal balance constraints
Q CHP , t + Q o u t , t + Q E , t = Q L , t + Q H , l o s s , t
where Q CHP , t represents the thermal output of the CHP unit at time t; Q L , t represents the thermal load of the system at time t; and Q H , l o s s , t represents the thermal loss of the system at time t.
(3) Wind power and photovoltaic unit output constraints
P WP , max P WP , t 0
P PV , max P PV , t 0
where P WP , max and P PV , max represent the maximum output of the WP and PV in period t, respectively.
(4) ETSB constraints
Through electric–thermal conversion and coupling, the ETSB can achieve peak cutting and valley filling of the load to reduce energy waste. Its operational state must meet the following constraints.
(1)
Capacity constraints
Based on the abovementioned ETSB model, the thermal storage state of the ETSB Q s t o r e , t should be less than its thermal storage capacity Q s t , c a at any time.
0 Q s t o r e , t Q s t , c a
(2)
Charge and discharge constraints
Restricted by planning solutions, the charge and discharge parameters of the ETSB must be kept within a certain range.
{ 0 P E , t P E , max 0 Q o u t , t Q o u t , max

3.2. Energy Flow Allocation of Lower Level

Obtaining the operating status of system equipment is a prerequisite for conducting economic analyses. Therefore, the energy distribution results of the system should be calculated at the lower level of optimal operation.

3.2.1. Objective Function

Taking the lowest operation cost of system as the objective function, we calculate the system energy flow when configuring the ETSB at different CHP nodes.
min C o p = C o p , p r + C o p , o m

3.2.2. Constraint Function

(1) Power grid constraints
Power grid constraints mainly include power flow constraints, voltage amplitude constraints, and line current constraints of the power grid, as follows:
  • Power flow constraints
In the operation of the power system, it is necessary to meet power flow constraints. The DistFlow power flow equation of distribution power networks is used in this study [20].
i v ( j ) [ P i j P i j 2 + Q i j 2 U i 2 R i j ] = l u ( j ) P j l + P j j Ω N
i v ( j ) [ Q i j P i j 2 + Q i j 2 U i 2 X i j ] = l u ( j ) Q j l + Q j j Ω N
U j 2 = U i 2 2 ( R i j P i j + X i j Q i j ) + [ R i j 2 + X i j 2 ] P i j 2 + Q i j 2 U i 2 i j Ω L
where u ( j ) and v ( j ) are the collection of nodes connected to node j, which are located downstream/upstream of node j; P i j and Q i j are the active power and reactive power flowing through branch i-j, respectively; U i and U j are the voltage amplitude of nodes i and j in their corresponding discrete states, respectively; P j l and Q j l are equivalent active and reactive loads at node j, respectively; R i j and X i j are the resistance and reactance of branch i-j, respectively; and Ω N and Ω L are the collection of nodes and branches of the system.
  • Voltage amplitude constraints
In the operational stage of the IES, the voltage of the power subsystem nodes needs to meet the following constraints.
U i min U i , t U i max
  • Power line current constraints
Due to the limitations of the line transmission capacity, the branch transmission current needs to satisfy certain constraints.
I i j 2 = P i j 2 + Q i j 2 U i 2 i j = Ω L
| I i j | I i j , max i j Ω L
  • CHP output constraints
P CHP , i max P CHP , t P CHP , i min
where P CHP , i max and P CHP , i min are the upper and lower limits of the power output of the CHP unit, respectively.
(2) Heating network constraints
Based on the quality regulation mode, the following constraints need to be met during the energy flow calculation process to ensure the safe and stable operation of the heating network [21].
  • Nodal flow rate balance constraints
According to Kirchhoff’s law, the flow rate of hot water entering and exiting nodes in the supply and return water network must be equal.
j S j p i p e + q j i n = j S j p i p e q j o u t
where q j i n and q j o u t are the water supply and return flow rate of node j, respectively, while S j p i p e + and S j p i p e are the pipe collection of the input and output node j, respectively.
  • Nodal thermal balance constraints
The thermal release energy of a node or line is directly proportional to the temperature and flow rate of hot water passing through it.
Q j = C q j ( T j i n T j o u t )
where Q j is the thermal release energy of node j; C is the specific heat capacity of water; q j is the water flow rate of node j; and T j i n is the water temperature of the input and output node j.
  • Pipeline temperature drop constraints
T i , e n d = ( T i , s t a r t T a ) e λ L i C q i
where T i , s t a r t and T i , e n d are the input and output temperature of pipeline i, respectively; T a is the ambient temperature; q i is the water flow rate of pipeline i; λ is the total heat transfer coefficient per unit length of pipeline; and L i is the length of pipeline i.
  • Nodal temperature constraints
According to the heating requirements, the temperature of the node’s supply and return water needs to meet the temperature constraints.
T i n , min T i i n T i n , max
T o u t , min T i o u t T o u t , max
where T i i n and T i o u t are the supply and return water temperature of node i, respectively; T i n , min and T i n , max are the supply water temperature upper and lower limit of node i, respectively; and T i n , min and T i n , max are the return water temperature lower and upper limit of node i, respectively.
  • Intersection node temperature constraints
Water flows through multiple pipelines and enters node k, which can be named the hydraulic intersection node. The temperature of the water flowing out of node k should satisfy the following constraints [22].
q k o u t T k o u t = n S k p i p e + ( q n i n T n i n )
where q n i n is the mass flow rate of pipeline n into the hydraulic intersection node k; q k o u t is the output mass flow rate of node k; T n i n is the input fluid temperature of pipeline n to node k; T k o u t is the output temperature of node k; S k p i p e + is the pipe collection of the input node k.

3.3. Model Solving Method

In view of the proposed bi-level planning and configuration model of ETSB, the upper level is a typical mixed-integer linear programming (MILP) problem, while the lower-level energy flow calculation has nonlinear characteristics. Thus, MILP and the genetic algorithm (GA) are combined in this study to solve the model [23].

3.3.1. Solution of the Upper Level Based on MILP

The goal of upper-level planning is to obtain the lowest energy supply cost for the system for different configuration nodes and parameters by solving the MILP process. The MILP problem presented in this paper can be described as follows:
min C T x s . t . { A x b A e q x = b e q x min x i x max x j { 0 , 1 }
where C T x is the objective function; A and Aeq are coefficient matrices; b and beq are n-dimensional coefficient vectors; and x max and x min are the upper and lower limits of x, respectively.
MILP combines the characteristics of linear and integer programming and has been widely applied in comprehensive energy system planning. Based on MILP, the upper-level planning process is as follows:
(1)
Step 1: Input the device parameters, load data, WP–PV power, and fuel price as initial data into the optimization model.
(2)
Step 2: Define the output of all the devices of EH-IES as decision variables.
(3)
Step 3: Set the objective function based on Equations (3)–(8) and set constraints based on results transferred from the lower level and Equations (9)–(14).
(4)
Step 4: Use the Cplex solver to solve the above optimization problem and obtain the configuration result.

3.3.2. Solution of the Lower Level Based on the GA

The characteristics of the multi-energy coupling of IES make its energy flow equation high-dimensional and strongly nonlinear. This also means that optimization problems containing multi-energy flow equations are essentially nonconvex and nonlinear optimization problems, which are difficult and time-consuming to solve directly. The GA has a good global search performance in solving nonlinear problems, so it is used to solve the optimal energy flow problem in the lower level. The calculation process is shown in Figure 4.
Based on the GA, the process of lower-level optimization is as follows:
(1)
Step 1: Input and initialize the system parameters based on the original parameters and the configuring results of the upper level.
(2)
Step 2: Generate the initial population and set the corresponding population size. In this paper, the initial population size is 400, the number of iteration steps is 100, the crossover rate is 0.85, and the mutation rate is 0.1.
(3)
Step 3: Evaluate the fitness of the population and generate a new population through selection, crossover, and mutation.
(4)
Step 4: Select two individuals from the parent generation for crossover mutation and add the generated new individuals to the population.
(5)
Step 5: Obtain the optimal solution by selecting the individual with the highest fitness during the calculation process and output the result.

4. Case Study

Taking an EH-IES at an industrial park in Ordos, China, as an example, the effectiveness of the proposed method was evaluated through analyzing the operational status of the EH-IES after configuring the ETSB. The original EH-IES without an ETSB is composed of both a nine-node electrical subsystem and an eight-node heating subsystem, containing two CHP units and one WP–PV joint output node, whose topology is shown in Figure 5.
It should be noted that the two CHP units in this system are a gas turbine unit (CHP 1) and an extraction condensing unit (CHP 2). Therefore, the relationship between the electrical and thermal energy of the two CHP units can be described as follows.
(1)
Gas turbine CHP unit
The relationship between electrical and thermal energy generation of the gas turbine CHP unit can be simplified into the following equation. Evidently, this indicates that the thermal–electrical energy ratio is constant, and it is set to 1.3 in this paper.
λ C H P , g = Q C H P , g P C H P , g
where λ C H P , g is the thermal–electrical energy ratio of the gas turbine CHP unit, while Q C H P , g and P C H P , g are the quantity of thermal and electrical energy generated by the gas turbine CHP unit, respectively.
(2)
Extraction condensing CHP unit
The relationship between electrical and thermal energy generation of the extraction condensing CHP unit is represented using a feasible region [24]. The area enclosed by the solid line and the Y-axis of Figure 6 represents the feasible operating region of the extraction condensing CHP unit used in this case.
It is clear that the extraction condensing CHP unit has good thermal–electrical regulation performance and can be combined with the gas turbine CHP unit to meet the energy supply requirements of the EH-IES in island mode.

4.1. Basic Data

Table 1 lists the resistance and reactance of the power subsystem lines, providing basic parameters for subsequent calculations and the analysis of power flow.
Table 2 lists the pipe length and diameter of the heat subsystem, providing basic parameters for subsequent calculations and the analysis of thermal energy flow.
Due to the configuration of the ETSB being optimized in an existing EH-IES, the parameters of the devices in the system (which are listed in Table 3) will have a significant impact on the calculation results.
According to the above mathematical mode, there are two parts to the ETSB, namely the EB and heat energy storage devices, and the basic parameters of the two parts are listed in Table 4.
Typical daily electrical and thermal loads are the foundation for optimizing resource configuration and improving overall energy supply efficiency. They can also help us to better understand the peak and off-peak load periods of the system. Due to the fact that electric and thermal loads are mainly used throughout the year in the industrial park, the whole year can be divided into the heating period (from 15 October to 15 April) and the non-heating period (from 16 April to 14 October) to more accurately evaluate the effect of ETSB in the system. The typical daily electrical and thermal loads of the two periods are determined by calculating the average value of each time during the corresponding period, as shown in Figure 7.
As seen in Figure 7, there is a significant increase in electrical load during the non-heating period compared to the heating period, and the peak–valley difference is not clear. In the heating period, the thermal load is typically high at night and low during the day, while the industrial heat demand is relatively stable in the non-heating period.
Additionally, in the planning and analysis process, the typical power output of the WP and PV is calculated from the average wind speed and solar irradiance at each time during the corresponding period, as shown in Figure 8.

4.2. Optimal Configuration Results

To expedite calculation of in this study, Matlab (R2021a) software suites are utilized to write computational programs based on the above basic data, and the analysis of results is as follows.

ETSB Parameter Results

The parameters of the maximum input power P E , t , the maximum discharge heat Q o u t , t , and the thermal storage capacity of ETSB Q s t , c a can not only greatly influence the energy flow allocation of the EH-IES but also influence the initial investment cost of the ETSB. When the configuration location of the ETSB is determined, the planning results will display significant differences due to the differences in these three parameters.
To provide a more detailed description of the changes, the double-layer optimization model was programmed using MATLAB software and the case solution was completed. Figure 9 shows the changes in optimal energy supply cost over a whole year, calculated by the planning model when the ETSB is configured at the CHP 2 node. Moreover, due to the specification parameters of the actual ETSB being designed in increments of 0.5 MW per level, the parameters of the ETSB are taken as 0.5 MW, 1.0 MW, 1.5 MW, 2.0 MW, and 2.5 MW.
When P E , max increases, it means that more surplus electricity can be converted to heat energy by the ETSB and participates in the heating supply, thereby reducing the system’s energy supply cost. But when P E , max exceeds a certain value, the constraints of investment and maintenance costs will be further amplified, resulting in an upward trend in the system’s energy supply cost. Therefore, the energy supply cost tends to increase after a decrease in any one curve in Figure 9, where Q s t , c a , max and Q o u t , max are fixed. The same reason can also be used to describe the changes in the energy supply cost with Q s t , c a , max . So, for five nodes of fixed P E , max values in any sub-graph, the energy supply cost also displays the tendency to decline at the beginning and rise later. Since there is only an impact on the thermal supply of the ETSB, the energy supply cost does not show a trend similar to the above two patterns with the increase in Q o u t , max . However, when P E , max and Q s t , c a , max are fixed, there still exists an optimal Q o u t , max value to make the heating supply more reasonable.
When configuring the ETSB at the CHP1 node, the same regulation can be obtained through the above calculation and analysis process, which is not described again here. Finally, the lowest energy supply cost for a whole year can be obtained by comparing the energy supply costs under each condition, and the corresponding configuration result is the optimal scheme for the ETSB, which is shown in Table 5.

4.3. Energy Balance Analysis

To facilitate further analyses and discussions of the effect of the ETSB on the EH-IES, the balance in energy supply and demand of the system during a typical day in the non-heating period and the heating period is shown in Figure 10 and Figure 11, in which the energy results of the original scheme (without an ETSB) are also presented for comparison purposes.

4.3.1. The Power Balance

In Figure 10a,b, the power balance during the non-heating and heating periods in the original scheme is presented, while Figure 10c,d present the power balance during the non-heating and heating periods in the optimal scheme, respectively.
In the non-heating period, the power output of both schemes displays identical values and can meet the load demands. The only difference is that the ETSB absorbs the excess electricity from the 23rd and 24th hours.
In the heating period, the power output of both schemes can also meet the load demands. However, the electrical and thermal loads do not match those shown in Figure 6, meaning that the thermal load is much higher in the heating period. Therefore, when the thermal load is high in the nighttime, the electrical output also increases accordingly, and a large amount of surplus electricity is produced. Due to the insufficient flexibility of the system, the surplus electricity can only be abandoned in the original scheme. Moreover, a clear reduction in the generation of electricity in the original scheme can be observed, especially during the daytime. More generated electricity means more surplus electricity. However, based on the electrical–thermal energy conversion characteristics of the ETSB in the optimal scheme, the surplus electrical energy can be converted into thermal energy to power heating or stored in the TSS during the low-electrical-load period. Therefore, the calculation results indicate that about 55.9 MWh of electrical energy is generated during a typical day in the heating period in the original scheme, while in the optimal scheme, the generated quantity increases to about 71.6 MWh, with about 26.8% or 19.6 MWh being consumed by the ETSB, effectively increasing the flexibility of system regulation.

4.3.2. The Thermal Balance

In Figure 11a,b, the thermal balance during the non-heating and heating periods in the original scheme is presented, while Figure 11c,d present the power balance during the non-heating and heating periods in the optimal scheme, respectively.
In the non-heating period, the total thermal generation of the two CHP units reaches 51.4 MWh, with 35.0 MWh being generated by CHP 1 and 16.4 MWh being generated by CHP 2. Meanwhile, the total thermal energy generation of the system drops to 45.4 MWh in the optimal scheme, with 38.5 MWh being generated by CHP 1 and 6.4 MWh being generated by CHP 2. Although the ETSB serves as a heat source, the amount of surplus energy is relatively small, so its utilization rate is not high. The thermal output of the ETSB is only 2.0 MW.
In the heating period, the total thermal generation of the original scheme increases to 81.7 MWh to satisfy the load, with 32.5 MWh generated by CHP 1 and 49.2 MWh generated by CHP 2. In addition to heating by two CHP units, the ETSB is also added as a heat source in the optimal scheme. A total of 79.4 MWh of thermal energy is generated during a typical day, with 37.9 MWh being generated by CHP 1 and 24.3 MWh being generated by CHP 2. Additionally, the ETSB plays an important role in the thermal supply process of the system as the third heat source, with the thermal output being 17.1 MWh, accounting for 21.5% of the total thermal generation, which reflects the energy spatiotemporal transfer characteristics of the ETSB.
The phenomenon of thermal abandonment only occurs in the original scheme, and the extent is low. Furthermore, it can be noted that the amount of surplus thermal energy is much less than the amount of surplus electricity, as shown in Figure 11. Combining this result with the fact that thermal load is much higher than the electrical load, as seen in Figure 6, the power abandonment phenomenon is not flexible, while the thermal energy demand is met, leading to a higher electricity output.
Additionally, due to the low thermal load during the non-heating season, there is a surplus of thermal energy generated by the system at the 2nd, 3rd, 14th, 15th, and 22nd hour. Meanwhile, because of the higher thermal load in the heating period, the surplus energy is all electricity. Thus, the phenomenon of surplus thermal energy being stored only appears in the results in Figure 11c.

4.4. Abandoned Energy Statistics and Analysis

The reason for energy abandonment is the strong coupling between electrical and thermal energy, which results in insufficient flexibility in system regulation. The ETSB can alleviate and improve the problem through electrical to thermal energy conversion technology and heat energy storage technology. To further illustrate the role of the ETSB in reducing power abandonment, Figure 12 shows the quantity of energy abandoned in the original scheme.
As shown in Figure 12, approximately 2.3 MWh of energy is abandoned during a typical day in the non-heating period, of which the majority is thermal energy, accounting for 92.2% of the total abandoned energy. Meanwhile, in the heating period, a total of 14.9 MWh of energy is abandoned within a day, comprising 5.6 MWh of thermal energy and 9.3 MWh of electricity. Moreover, it is obvious that the phenomenon of electricity abandonment occurs at the 10th, 11th, and 16–19th hours, corresponding to the low-electrical-load periods within a day.
Electric–thermal decoupling is achieved and energy is spatiotemporally transferred by configuring the ETSB in the system. Therefore, due to the improved adjustment flexibility of the system after configuring the ETSB, the phenomenon of energy abandonment is eliminated.

4.5. Energy Storage of ETSB

In addition to the ability to convert electrical energy to thermal energy, ETSB also has the function of storing thermal energy. Figure 13 shows the thermal storage state of the ETSB in the optimal scheme.
In Figure 13, it can be seen that the ETSB plays a significant role in heat storage both during the heating and non-heating periods. In the non-heating period, there are three obvious peaks in the amount of thermal energy stored in the ETSB, appearing at the 3rd, 15th, and 22nd hour, with the thermal energy capacity of 0.9 MWh and 0.8 MWh, corresponding to the time at which surplus thermal energy is present. In the heating period, due to the large amount of surplus electricity, the ETSB’s participation in electric–thermal energy conversion and heating becomes more frequent. Therefore, the thermal storage states are maintained at a high level throughout the day. The highest thermal storage capacity reaches 1.6 MWh at the 15th hour.

4.6. The Energy Flow Allocation Analysis

One of the methods to implement the above functions of the ETSB is to change the energy flow allocation of the system. Contrastive analysis of the variation in energy flow allocation after configuring the ETSB in EH-IES is helpful to understand this effect. Figure 14 shows the energy flow allocation of the energy lines of the system at the 8th hour of the heating period in the optimal and original schemes [25].
As shown in Figure 14, about 0.64 MW of electricity generated by CHP 2 is abandoned in the original scheme because it is the balance node of the system. But after configuring the ETSB in the system, it consumes about 0.45 MW of surplus electricity in this hour, which causes the electrical output of CHP 1 and CHP 2 to change from 1.32 MW and 0.41 MW to 1.61 MW and 0.58 MW, respectively. Subsequently, the electrical allocation of the entire subsystem’s lines varies, not only referring to the values, but also indicating the direction of the power flow, such as the power line between Nodes 2 and 8.
On the other hand, the thermal output of CHP 1 and CHP 2 also changes from 1.72 MW and 1.89 MW to 2.09 MW and 1.05 MW, respectively, after configuring the ETSB due to the coupling of electricity and heat of CHP units. Correspondingly, the thermal energy flow allocations of each pipeline change, with the pipeline between nodes 1 and 2 being particularly representative, with its thermal energy flow value changing from 0.08 MW to 0.3 MW and its direction being reversed.

5. Conclusions and Future Work

In this paper, an ETSB is introduced into an EH-IES to further enhance system flexibility and ensure the reliability of the energy supply, and a bi-level optimal configuration model based on the established model of ETSB is proposed. An example EH-IES is simulated to obtain the optimal planning and configuration of ETSB. By comparing the operation results of the original scheme and the scheme after configuring the ETSB, the impact of the ETSB is analyzed, and the effectiveness of the proposed planning scheme is verified. According to the simulation results, the following conclusions are obtained:
(1)
The energy supply cost of the EH-IES is closely linked to the selection of the ETSB parameters, namely P E , t , Q o u t , t , and Q s t , c a .
(2)
The lack of flexibility in the original scheme is the main reason for energy abandonment. After the ETSB is optimally configured, the flexibility of system regulation is greatly increased and the phenomenon of energy abandonment is eliminated.
(3)
The main method through which ETSB improves system flexibility is to redistribute energy within the system by increasing the electric–thermal decoupling ability and energy spatiotemporal transfer, thus optimizing the electrical/thermal source output and energy flow distribution.
Research into the optimal configuration of ETSB in EH-IES was conducted, and the relevant system energy flow allocation has been included in the analysis. However, the data description process is at a relatively early stage and does not consider the amount of carbon emissions. In future work, we plan to introduce energy and carbon flow tracking methods to thoroughly analyze the impact of ETSB on IES energy and carbon flows to further analyze their role in reducing carbon emissions.

Author Contributions

Conceptualization, X.Z. (Xiaoming Zhang) and P.Y.; methodology, X.Z. (Xiaoming Zhang); software, J.F. and G.L.; formal analysis, X.Z. (Xiaoming Zhang) and J.F.; investigation, C.D.; writing—original draft preparation, J.F.; writing—review and editing, X.Z. (Xiaoming Zhang); supervision, X.Z. (Xin Zhang); funding acquisition, X.Z. (Xiaoming Zhang). All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Science and Technology Planning Project of the Inner Mongolia Autonomous Region under Grant No. 2022YFHH0027 and 2023YFHH0055, Inner Mongolia Natural Science Foundation under Grant No. 2022LHQN05002, Grassland Talent Project High-level training Project of the Inner Mongolia Autonomous Region under Grant No. CYYC012054, the central government guides local science and technology development fund project under Grant No. 2023ZY0020, and basic research funds for universities directly under the Inner Mongolia Autonomous Region under Grant No. 2023QNJS199.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Zhang, D.; Tian, J.; Goh, H.-H.; Liu, H.; Li, X.; Zhu, H.; Wu, X. The Key Technology of Smart Energy System and Its Disciplinary Teaching Reform Measures. Sustainability 2022, 14, 14207. [Google Scholar] [CrossRef]
  2. Li, J.; Xu, W.; Cui, P.; Qiao, B.; Feng, X.; Xue, H.; Wang, X.; Xiao, L. Optimization Configuration of Regional Integrated Energy System Based on Standard Module. Energy Build. 2020, 229, 110485. [Google Scholar] [CrossRef]
  3. Lu, S.; Li, Y.; Xia, H. Study on the configuration and operation optimization of CCHP coupling multiple energy system. Energy Convers. Manag. 2018, 177, 773–791. [Google Scholar] [CrossRef]
  4. Huang, W.; Zhang, N.; Yang, J.; Wang, Y.; Kang, C. Optimal configuration planning of multi-energy systems considering distributed renewable energy. IEEE Trans. Smart Grid 2019, 10, 1452–1464. [Google Scholar] [CrossRef]
  5. Qiao, Y.; Hu, F.; Xiong, W.; Guo, Z.; Zhou, X.; Li, Y. Multi-objective optimization of integrated energy system considering installation configuration. Energy 2023, 263, 125785. [Google Scholar] [CrossRef]
  6. Meibom, P.; Kiviluoma, J.; Barth, R.; Brand, H.; Weber, C.; Larsen, H.V. Value of electric heat boilers and heat pumps for wind power integration. Wind. Energy 2007, 10, 321–337. [Google Scholar] [CrossRef]
  7. Zhao, P.; Gou, F.; Xu, W.; Shi, H.; Wang, J. Energy, exergy, economic and environmental (4E) analyses of an integrated system based on CH-CAES and electrical boiler for wind power penetration and CHP unit heat-power decoupling in wind enrichment region. Energy 2023, 263, 125917. [Google Scholar] [CrossRef]
  8. Wang, W.; Yuan, B.; Sun, Q.; Wennersten, R. Application of energy storage in integrated energy systems—A solution to fluctuation and uncertainty of renewable energy. J. Energy Storage 2022, 52, 104812. [Google Scholar] [CrossRef]
  9. Bazdar, E.; Sameti, M.; Nasiri, F.; Haghighat, F. Compressed air energy storage in integrated energy systems: A review. Renew. Sustain. Energy Rev. 2022, 167, 112701. [Google Scholar] [CrossRef]
  10. Ujjwal, D.; Akhtar, K.; Juan, S. A review of key functionalities of battery energy storage system in renewable energy integrated power systems. Energy Storage 2021, 3, 224. [Google Scholar]
  11. Ren, H.; Jiang, Z.; Wu, Q.; Li, Q.; Yang, Y. Integrated optimization of a regional integrated energy system with thermal energy storage considering both resilience and reliability. Energy 2022, 261, 125333. [Google Scholar] [CrossRef]
  12. Wei, W.; Guo, Y.; Hou, K.; Yuan, K.; Song, Y.; Jia, H.; Sun, C. Distributed Thermal Energy Storage Configuration of an Urban Electric and Heat Integrated Energy System Considering Medium Temperature Characteristics. Energies 2021, 14, 2924. [Google Scholar] [CrossRef]
  13. Böttger, D.; Götz, M.; Theofilidi, M.; Bruckner, T. Control power provision with power-to-heat plants in systems with high shares of renewable energy sources—An illustrative analysis for Germany based on the use of electric boilers in district heating grids. Energy 2015, 82, 157–167. [Google Scholar] [CrossRef]
  14. Sinha, R.; Jensen, B.; Pillai, J. Impact Assessment of Electric Boilers in Low Voltage Distribution Network. In Proceedings of the 2018 IEEE Power & Energy Society General Meeting, Portland, OR, USA, 5–10 August 2018; pp. 1–5. [Google Scholar]
  15. Zhao, S.; Ge, Z.; Sun, J.; Ding, Y.; Yang, Y. Comparative study of flexibility enhancement technologies for the coal-fired combined heat and power plant. Energy Convers. Manag. 2019, 184, 15–23. [Google Scholar] [CrossRef]
  16. Lei, Z.; Wang, G.; Li, T.; Cheng, S.; Yang, J.; Cui, J. Strategy analysis about the active curtailed wind accommodation of heat storage electric boiler heating. Energy Rep. 2021, 7, 65–71. [Google Scholar] [CrossRef]
  17. Liu, B.; Li, J.; Zhang, S.; Gao, M.; Ma, H.; Li, G.; Gu, C. Economic Dispatch of Combined Heat and Power Energy Systems Using Electric Boiler to Accommodate Wind Power. IEEE Access 2020, 8, 41288–41297. [Google Scholar] [CrossRef]
  18. Yang, L.; Zhang, X.; Gao, P. Research on heat and electricity coordinated dispatch model for better integration of wind power based on electric boiler with thermal storage. IET Gener. Transm. Distrib. 2018, 12, 3736–3743. [Google Scholar] [CrossRef]
  19. Li, P.; Wang, Z.; Liu, H.; Wang, J.; Guo, T.; Yin, Y. Bi-level optimal configuration strategy of community integrated energy system with coordinated planning and operation. Energy 2021, 236, 121539. [Google Scholar] [CrossRef]
  20. Baran, M.; Wu, F. Network reconfiguration in distribution systems for loss reduction and load balancing. IEEE Trans. Power Deliv. 1989, 4, 1401–1407. [Google Scholar] [CrossRef]
  21. Ma, H.; Liu, C.; Zhao, H.; Zhang, H.; Wang, M.; Wang, X. A novel analytical unified energy flow calculation method for integrated energy systems based on holomorphic embedding. Appl. Energy 2023, 344, 121163. [Google Scholar] [CrossRef]
  22. Hassine, I.; Eicker, U. Impact of load structure variation and solar thermal energy integration on an existing district heating network. Appl. Therm. Eng. 2013, 50, 1437–1446. [Google Scholar] [CrossRef]
  23. Khamees, A.; Badra, N.; Abdelaziz, A. Optimal power flow methods: A comprehensive survey. Int. Electr. Eng. J. 2016, 7, 2228–2239. [Google Scholar]
  24. Sondergren, C.; Ravn, H.F. A method to perform probabilistic production simulation involving combined heat and power units. IEEE Trans. Power Syst. 1996, 11, 1031–1036. [Google Scholar] [CrossRef]
  25. Cheng, Y.; Zhang, N.; Wang, Y.; Yang, J.; Kang, C.; Xia, Q. Modeling Carbon Emission Flow in Multiple Energy Systems. IEEE Trans. Smart Grid 2019, 10, 3562–3574. [Google Scholar] [CrossRef]
Figure 1. Framework diagram of an EH-IES.
Figure 1. Framework diagram of an EH-IES.
Electronics 13 03567 g001
Figure 2. Energy flow diagram of ETSB.
Figure 2. Energy flow diagram of ETSB.
Electronics 13 03567 g002
Figure 3. The expansion planning scheme of ETSB.
Figure 3. The expansion planning scheme of ETSB.
Electronics 13 03567 g003
Figure 4. Flowchart for the GA calculation of the lower level.
Figure 4. Flowchart for the GA calculation of the lower level.
Electronics 13 03567 g004
Figure 5. Topology of the original EH-IES.
Figure 5. Topology of the original EH-IES.
Electronics 13 03567 g005
Figure 6. The operational feasible region of the extraction condensing CHP unit.
Figure 6. The operational feasible region of the extraction condensing CHP unit.
Electronics 13 03567 g006
Figure 7. Electrical and thermal load demand curve. (a) Non-heating period; (b) heating period.
Figure 7. Electrical and thermal load demand curve. (a) Non-heating period; (b) heating period.
Electronics 13 03567 g007
Figure 8. Typical daily WP−PV output curve. (a) Non−heating period; (b) heating period.
Figure 8. Typical daily WP−PV output curve. (a) Non−heating period; (b) heating period.
Electronics 13 03567 g008
Figure 9. The relationship of optimal results and ETSB parameters. (a) Q o u t , max = 0.5 ; (b) Q o u t , max = 1.0 ; (c) Q o u t , max = 1.5 ; (d) Q o u t , max = 2.0 ; (e) Q o u t , max = 2.5 ; (f) Q o u t , max = 3.0 .
Figure 9. The relationship of optimal results and ETSB parameters. (a) Q o u t , max = 0.5 ; (b) Q o u t , max = 1.0 ; (c) Q o u t , max = 1.5 ; (d) Q o u t , max = 2.0 ; (e) Q o u t , max = 2.5 ; (f) Q o u t , max = 3.0 .
Electronics 13 03567 g009aElectronics 13 03567 g009b
Figure 10. Power balance of the two schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme; (c) Non−heating period of the optimal scheme; (d) Heating period of the optimal scheme.
Figure 10. Power balance of the two schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme; (c) Non−heating period of the optimal scheme; (d) Heating period of the optimal scheme.
Electronics 13 03567 g010
Figure 11. Thermal balance of the two schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme; (c) Non−heating period of the optimal scheme; (d) Heating period of the optimal scheme.
Figure 11. Thermal balance of the two schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme; (c) Non−heating period of the optimal scheme; (d) Heating period of the optimal scheme.
Electronics 13 03567 g011aElectronics 13 03567 g011b
Figure 12. The abandoned electricity in both schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme.
Figure 12. The abandoned electricity in both schemes. (a) Non−heating period of the original scheme; (b) Heating period of the original scheme.
Electronics 13 03567 g012
Figure 13. Thermal storage state of the ETSB in the optimal scheme.
Figure 13. Thermal storage state of the ETSB in the optimal scheme.
Electronics 13 03567 g013
Figure 14. The energy flow allocation of the system. (a) Original scheme; (b) Optimal scheme.
Figure 14. The energy flow allocation of the system. (a) Original scheme; (b) Optimal scheme.
Electronics 13 03567 g014
Table 1. The basic parameters of the electrical network.
Table 1. The basic parameters of the electrical network.
NumberInitial NodeFinal NodeR (ohm)X (ohm)
1160.0426400.0208
2120.0026240.0128
3280.0032800.0160
4380.0052480.0256
5340.0037720.0184
6470.0027880.0136
7570.0042640.0208
8690.0008000.0120
Table 2. The basic parameters of the heat pipeline.
Table 2. The basic parameters of the heat pipeline.
NumberInitial NodeFinal NodePipe Length (m)Pipe Diameter (mm)
181257.6125
21497.540
3125140
425595100
523271.332
636235.465
773177.340
Table 3. The basic parameters of the energy devices.
Table 3. The basic parameters of the energy devices.
EquipmentThe Bounds Value (MW)Power Generation EfficiencyMaintenance Cost(CNY/kW)
CHP10~2.50.30.025
CHP20~1.50.250.025
WP0~0.8-0.0196
PV0~0.8-0.0235
Table 4. The parameters of the ETSB.
Table 4. The parameters of the ETSB.
EquipmentEfficiencyInvestment
(CNY/kW)
Operation Cost
(CNY/kW)
Lifetime
(Year)
EB0.95500.0415
TSS0.93000.0215
Table 5. The optimal scheme of the case.
Table 5. The optimal scheme of the case.
Configuration LocationPE,t (MW)Qst,ca (MWh)Qout,t (MW)Energy Supply Cost (CNY)
CHP 21.52.01.019,568,068
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhang, X.; Feng, J.; Liang, G.; Ding, C.; Yang, P.; Zhang, X. Bi-Level Optimal Configuration of Electric Thermal Storage Boilers in Thermal–Electrical Integrated Energy System. Electronics 2024, 13, 3567. https://doi.org/10.3390/electronics13173567

AMA Style

Zhang X, Feng J, Liang G, Ding C, Yang P, Zhang X. Bi-Level Optimal Configuration of Electric Thermal Storage Boilers in Thermal–Electrical Integrated Energy System. Electronics. 2024; 13(17):3567. https://doi.org/10.3390/electronics13173567

Chicago/Turabian Style

Zhang, Xiaoming, Jiaoyang Feng, Guangzhe Liang, Chonglei Ding, Peihong Yang, and Xin Zhang. 2024. "Bi-Level Optimal Configuration of Electric Thermal Storage Boilers in Thermal–Electrical Integrated Energy System" Electronics 13, no. 17: 3567. https://doi.org/10.3390/electronics13173567

APA Style

Zhang, X., Feng, J., Liang, G., Ding, C., Yang, P., & Zhang, X. (2024). Bi-Level Optimal Configuration of Electric Thermal Storage Boilers in Thermal–Electrical Integrated Energy System. Electronics, 13(17), 3567. https://doi.org/10.3390/electronics13173567

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop