Optimal high‐level control of building HVAC system under variable price framework using partially linear model

National Key R&D Programme of China, Grant/ Award Number: 2020YFB0905900 Abstract As typical thermostatic loads, heating, ventilation, and air conditioning (HVAC) systems possess high flexibility due to buildings' thermal inertia. It is beneficial to coordinate HVAC loads with power grid operations. However, it is difficult to model the HVAC loads of realistic buildings. A partially linear model is proposed to capture the electrothermal characteristics of building HVAC loads. The proposed model combines the explainability of a physically based model and the fitting capability of a data‐driven model. The model retains linearity regarding mapping the control variables to the HVAC load and thus is compatible with existing studies based on the assumption of a succinct HVAC load model. Based on the proposed model, an optimal control strategy is developed to minimize the HVAC operational cost while considering comfort requirements of occupants. EnergyPlus is used as the building simulation engine to validate the proposed model and optimal control strategy.


| INTRODUCTION
With the fast process of global urbanization, the building sector has become a major energy consumer. Worldwide, building construction and operations together are responsible for nearly 40% of greenhouse gas emissions [1]. Among the appliances in a building, heating, ventilation, and air conditioning (HVAC) systems account for a large proportion of electricity consumption. According to a survey conducted by the U.S. Energy Information Administration, this proportion is approximately 48% [2].
In the field of building management, efforts have been paid to modelling and controlling HVAC loads in sophisticated ways so that energy saving and thermal comfort can both be improved. Simulation-based control methods search for the optimal strategy by running building energy simulations. In [3], a self-tuning load management strategy was designed to implement flexible control of the temperature setpoints in multiple zones. In [4], an energy management system for buildings in a microgrid was developed to optimize both energy consumption and thermal comfort, and the various occupancy schedules and fluctuations of solar energy were considered. Model predictive control (MPC) methods first establish the dynamic model of the system and then implement optimal control over a moving horizon. An energy-efficient MPC controller for HVAC operation was designed in [5] with regard to humidity and latent heat. In [6], a reliable HVAC load model and a hierarchical MPC structure were developed. In [7], recent studies on artificial neural network-based MPC were reviewed. Adaptive control methods enable flexible changes in the control gains so that the controllers can adapt to time-varying factors [8]. The advancement in adaptive control of building HVAC was reviewed in [9].
Proper management of building HVAC not only can reduce energy consumption but also can benefit power grids. As typical thermostatic loads, HVAC systems have the potential to provide flexibility to power grids due to buildings' thermal capacity. HVAC loads can be adjusted temporarily with a very limited impact on users' comfort. Large numbers of HVAC systems in urban building stock make this potential especially valuable for the economic operation of the power grid in urban areas. With appropriate demand-side management, optimization of HVAC system operation can greatly improve the energy efficiency of the grid with a high penetration of renewable energy. Therefore, modelling and optimal control of HVAC systems have aroused wide research interest.
Many studies have illustrated the application of HVAC systems in shedding/shifting loads, accommodating renewable energy, and providing ancillary services. For example, [10,11] performed optimal control of indoor temperature relative to real-time pricing. In [12], air-conditioning loads were modelled as virtual energy storage and were dispatched by the load aggregator to mitigate the volatility of wind power generation and hedge against economic loss in the real-time market. To mitigate fluctuations in PV generation, MPC strategies were proposed in [13] to coordinate a cluster of residential air conditioners. HVAC systems can also provide frequency and voltage regulation services to the power grid by following certain signals from the system operator [14,15].
However, the majority of existing literature have been based on the assumption that the electrothermal characteristics of building HVAC loads were known. The most commonly used method to describe HVAC loads is the equivalent thermal parameter (ETP) model. By analogy, the ETP model translates the thermal transfer process into an electric circuit diagram with elements representing thermal capacity and thermal resistances. Therefore, the ETP model is also referred to as the RC model. The complexity of the RC model depends on to what degree the inner heat transfer processes are ignored. The most simplified first-order RC model is composed of one capacitance representing the thermal capacity of the whole building and one resistance representing the thermal resistance of the outer wall, as applied in many studies. More sophisticated models use multiple capacitances and resistance components to depict the heat flows between different zones of the building.
Studies focussing on the interaction between building HVAC loads and the power grid usually assign some given values to ETP without considering how to determine the ETP values of real buildings, which has limited the practical application of these theories. Some studies have discussed monitoring and model parameter estimation of actual buildings. For example, in [16], dual estimation of unmeasured states and parameters of a building energy model was proposed based on Kalman filtering. In [17], a control performance monitoring method was implemented in a practical application. In [18], the heating demand of the building was estimated by combining autoregressive behaviour and longterm trends. The parameters of the HVAC load model can be estimated by analytical or experimental methods [19], but the former requires complete information of building structure, and the latter requires special experiments. In large-scale applications, both methods are intractable. For buildings in operation, HVAC loads are also largely affected by internal heat gains from occupants and electric appliances, as well as external heat gains via conduction, convection, and solar radiation.
To identify the electrothermal characteristics of real building HVAC loads, some studies have applied data-driven approaches for modelling and optimal control. In [20], several machine-learning techniques were compared for their performance of forecasting zonal cooling airflow rates, based on which HVAC loads were operated as a virtual battery. In [21], artificial neural networks were trained to forecast indoor temperatures of different zones in an office building, and optimal demand response of the HVAC system was implemented. Input convex recurrent neural networks were constructed in [22] to forecast HVAC loads, and the convexity significantly simplified building energy optimization problems. These data-driven methods have better forecast performance but are usually lacking in simplicity and physical explainability. Moreover, if we want to find an optimal control policy based on a black box model, the difficulty is greatly increased because of the lack of an explicit expression.
In summary, physically based RC models have been applied in many studies, but their parameters are hard to estimate because these parameters are influenced by complex and timevariant factors. The black box models may perform better to capture the realistic thermal dynamics of buildings, but these models are too complicated to be embedded in large-scale aggregation and optimization.
To solve this problem, a partially linearized HVAC model is proposed that combines the physically based model and a datadriven model. The partially linearized HVAC model consists of a linear and a non-linear part. The linear part is derived from the physical characteristics of HVAC loads, and the non-linear part applies data-driven methods to predict the complex effects of environmental factors. The proposed model captures the electrothermal characteristics of building HVAC loads from historical operating data. Additionally, it retains simplicity with regard to mapping the control variables to the electricity consumption of HVAC. Therefore, the proposed model is compatible with existing studies based on the first-order RC model. Based on the proposed model, an MPC strategy is developed to reach the optimal trade-off between the electricity bill and the occupants' thermal discomfort. The objective function is the weighted sum of the electricity cost and the thermal discomfort penalty, which is minimized by flexibly adjusting the indoor temperature setpoint under time-varying electricity prices.
The rest of this paper is organized as follows. Section 2 introduces the framework of the proposed model. Section 3 discusses the RC model and the HVAC system model to deduce the physical expression of HVAC loads. Section 3 illustrates the formulation of the partially linearized HVAC model, and an optimal control strategy is proposed in Section 4. The proposed model and strategy are validated by case studies in Section 5, and conclusions are discussed in Section 6.

| FRAMEWORK OF PROPOSED METHODS
The proposed methods aim to model building HVAC loads using a partially linear model and perform optimal control of the temperature setpoint considering dynamic electricity prices based on the load model. The structure of the proposed HVAC load model is shown in Figure 1. The left part is the grey box model derived from physical characteristics, and the right part is the black box model where data-driven methods are applied. Temperature setpoint, the controllable variable in the model, is the input of the grey box model. Therefore, the correlation between the temperature setpoint and the HVAC load is clear regardless of the complexity of the black box model. This correlation, represented by the parameters of the grey box model, is then used for optimal control.
The flow diagram of the proposed method is shown in Figure. 2 and consists of two main stages, the offline fitting stage and the online optimizing stage. In the offline fitting stage, historical operational data of the building are used to train an HVAC load model. Then, in the online optimizing stage, the load model is used to calculate the optimal temperature setpoint. The general framework of the proposed control method is identical with that of MPC, first establishing the dynamic model of the system and then implementing optimal control over a moving horizon. Slightly different from the typical MPC, the proposed control method is used to flexibly control the setpoint under a variable electricity price instead of tracking a desired output, and the proposed method consists of a data-driven component. In a wide sense, the proposed control method belongs to the domain of MPC.

| Thermal dynamics of buildings
As pointed out in [23], load modelling and short-term forecasting are essential for the implementation of automated load optimization. To perform optimal control of building HVAC loads while satisfying user comfort, the operator needs to know a building's electrothermal characteristics, specifically mapping the indoor temperature to the electricity consumption of HVAC. However, the realistic electrothermal characteristics of buildings are very complicated and hard to estimate. The most accurate method is to construct differential equation models of the building and its HVAC systems based on physical characteristics and then run simulations. Such white box models require complete information about the building and its HVAC systems and are typically solved with simulation software, such as TRNSYS [24] and EnergyPlus [25]. However, due to the lack of an explicit expression of electrothermal characteristics, large numbers of simulation runs are needed to search for the optimal control policy. This process is both time consuming and computationally expensive.
To simplify building electrothermal characteristics, the RC model is applied to describe the heat conduction process of building envelopes. Figure 3 depicts the basic first-order RC model of a thermal zone, where C represents the thermal capacity of the indoor air and furniture. R is the equivalent heat resistance of envelopes. T in is the indoor temperature. T out is the outdoor temperature. P cond is the heat conduction power flow caused by the temperature difference, and P T is the total effective heating/cooling rate. Hence, the differential equation describing the building's electrothermal characteristics is shown in Equation (1). Because the outdoor temperature changes much more slowly, it is treated as a constant.
When the building consists of multiple thermal zones with different temperatures, the first-order RC model can be extended to a network. The differential equations describing the temperature of each thermal zone are given by Equation (2), where the superscript i denotes the ith zone. A i is the set of zones adjacent to the ith zone. R i is the heat resistance between the ith zone and the outside environment, and R ij is the heat resistance between the ith and jth zones. The validity  of the RC network was tested in [26]. It should be noted that the total effective heating/cooling rate P T is the sum of the heating/cooling rate provided by the HVAC system and the heat gains from other means, including occupants, electric appliances, and solar irradiation:

| HVAC loads
The RC model only describes the heat conduction of the building envelopes. However, HVAC loads depend not only on the building's thermal characteristics but also on the characteristics of the HVAC system. The basic structure of an HVAC system serving N zones is shown in Figure 4. Assume q i ðtÞ is the airflow rate provided for the ith zone. The return air from all the zones is joined together, and the total airflow is denoted as q total ðtÞ. In the mixing box, part of the return air is replaced by fresh outdoor air. The temperature of the mixed air is shown in Equation (4), where δ is the recirculation ratio. The supply air passes through the chiller/heater, and its temperature changes to the setpoint T c , and then is sent to each zone by a fan. The electric power of the chiller/heater is shown in Equation (5), where the superscript E in P E c=h ðtÞ means electric power. c p is the heat capacity ratio of the air. COP is the coefficient of performance of the chiller/heater. Then, the supply air is compressed by a fan and sent to each zone. To keep the zonal temperature at the setting level, there is usually a thermostat in each zone that controls the supply airflow rate through a damper. The heating/cooling rate provided by the supply air for each zone can be calculated by Equation (6). The electric power of the fan is usually a cubic polynomial of the total airflow, as shown in Equation (7). The electric power of the entire HVAC system P E HVAC is the sum of P E c=h and P f an : P fan ðtÞ ¼ c 3 q 3 total ðtÞ þ c 2 q 2 total ðtÞ þ c 1 q total ðtÞ þ c 0 ð7Þ 4 | PARTIALLY LINEARIZED HVAC MODEL

| Preliminaries of partially linear model
In statistics, the additive model was first introduced by Friedman and Stuetzle in 1981 [27] and then extended to the generalized additive model by Hastie and Tibshirani in 1990 [28]. The basic structure of an additive model is shown in Equation (8), where the expectation of response variable EðY Þ is expressed as the sum of the effects of some predictors. The function f i can take any form, and x i can be either a vector or a scalar, which gives an additive model great adaptability and makes it less affected by the curse of dimensionality compared with a single high-dimensional function taking all the predictors as inputs.
If the additive model consists of both parametric components and non-parametric components, it is referred to as a semi-parametric model. A parametric model makes assumptions about the probability distribution of the data based on prior knowledge, and the model can be specified by determining a set of parameters that usually have physical meanings. In contrast, a non-parametric model requires no information about data distribution, and the functional form is determined F I G U R E 4 The basic structure of an HVAC system F I G U R E 3 First-order RC model by data using smoothers. Combining the advantages of both parametric and non-parametric models, semi-parametric models are especially suitable for situations where the effects of some predictors are known while the rest are not. In such situations, parametric models may have poor performance due to misspecified assumptions, whereas non-parametric models fail to take advantage of prior knowledge about the relationship between certain predictors and the response variable. A partially linear model is a commonly used form of semiparametric models, in which the response variable is supposed to have a linear relationship with a subset of predictors. It was first proposed by Engle et al. to study the relation between weather and monthly electricity sales [29]. The basic structure of a partially linear model is shown in Equation (9), where i ∈ f1; 2; …; ng is the order of the sample point, Y i is the vector of the remaining predictors, β ¼ ðβ 1 ; …; β p Þ T is a vector of linear coefficients, g is an unknown function to be estimated, and ε i is the error term.

End for End while
It is worth noting that the selection of the function to estimate g is quite flexible and should be based on the problem situation. For example, the estimator can be expressed as the sum of several smoothing functions or non-linear functions.
In a partially linear model, the linear term retains the simplicity and explainability, while the non-linear term accounts for the complex effects of other predictors. Therefore, the partially linear model strikes a balance between the accuracy and the usability. Because of this major advantage, partially linear models have been widely applied in statistics and economics. In these applications, researchers are usually interested in the contributions of particular explanatory variables, whereas other explanatory variables also impose their effects in different forms. For example, the combined effects of economic growth and urbanization on CO 2 emissions were investigated using a partially linear additive model in [30]. In [31], the impacts of renewable energy technology innovations on China's green productivity were examined.
In a partially linear model, if the estimator of the unknown function g is determined to be a certain form of smoother, the linear coefficient vector β and the function g can be estimated rather easily [32]. For example, define an estimator of g by Equation (10), where n is the total number of sample points, w nj ðt Þ are positive weight functions depending on the distance between t and T j . Assume β is known, and then gðT j Þ can be replaced by Y i − X T i β, as shown in Equation (11). Combining Equations (9) and (11), we obtain Equation (12). By rewriting is actually a linear regression model given the weight function w nj . Therefore, β can be obtained by applying the least square method, which is denoted as β LS . Replacing β by β LS , we obtain the estimator g n ðt Þ. When the estimator takes the form of other smoothers such as splines, similar processes can be performed: However, if a more complex function is used to approximate g, such as a decision tree or an additive model consisting of several smoothers, it cannot be expressed by β in the same way as in Equation (11). Under such circumstances, a more general method, namely a backfitting algorithm [28], should be applied. For an additive model defined by Equation (8), the backfitting algorithm, which is a heuristic approach, iteratively estimates each component.

| Partially linearized HVAC model
For a building served by one HVAC system, the temperature differences between zones are usually rather small. Therefore, the heat flows inside the building are ignorable compared with the heat flows from outside and the internal heat gains. We assume that the temperature of the whole building can be represented by one value, T in , so that the RC model and the HVAC model can be greatly simplified. Without loss of generalization, when the temperature differences in a building cannot be ignored, each thermal zone is taken as a 'building', and the adjacent zones are taken as the 'outside environment'.
By writing Equation (1) into time-discrete form, we have Equation (13), where the subscript k is the time-period index. The electric power of the fan is relatively small, so we simply ignore the non-linearity in Equation (7). Therefore, we can obtain from Equations (3)-(7) that the electric power of the HVAC system is in a linear relationship with the total heating/ cooling rate providing for the building, as shown in Equation (14). The total effective heating/cooling rate of the HE ET AL.
-217 building P T k is the sum of the supply heating/cooling rate P supply k and the heat gains from other ways P G k , as shown in Equation (15). Combining Equations (13)- (15), the electric power of the HVAC system is given by Equation (16). P G k is relevant to many factors, such as the density of occupants and solar irradiation, so it is hard to find an explicit expression for it. To solve this problem, we formulate it as an unknown nonparametric function of relevant factors and estimate the function by data-driven methods. Therefore, Equation (16) can be simplified into Equation (17), where d k is a vector of relevant environmental factors and f is the unknown function. This is a partially linear model and can be solved by methods mentioned in Section 4.1: To implement the HVAC load model, the relevant elements in the vector d k are chosen to be the solar irradiation rate, the humidity, the outside temperature, and the timestamp (current hour). The solar irradiation rate directly affects the building's heat gains, and hence we assume that the relationship between solar irradiation and HVAC load is linear. Humidity and outside temperature are related to the heat convection process. It is worth noting that although the outside temperature appears in the linear part, its effect is actually non-linear. There, the outside temperature is taken as the input of the nonparametric again. The occupant density and the working condition of electric appliances are related to time, so a timestamp is considered. The function f chosen to estimate the effects of the humidity, the outside temperature, and the timestamp is Gaussian process regression, which outperforms other commonly used regression methods in this study.
When the temperature setpoint of the building HVAC system is not optimized, the setpoint is usually kept at a baseline value T base that is the most suitable for the occupants. It should be noted that appropriate data are essential to training the HVAC load model. If the indoor temperature setting is kept constant during simulation, the parameters of the load model would be wrongly estimated. So we assume that the building has already adopted a simple stepwise strategy to change its temperature setpoint regarding the real-time price. The stepwise temperature setting strategy is illustrated by an example in Equation (18):

| OPTIMAL CONTROL STRATEGY
After obtaining the HVAC load model of the building, the temperature setpoint can be controlled to reach the optimal trade-off between the electricity bill and the occupants' thermal discomfort. When the setpoint is kept T base , the baseline load is calculated by Equation (19). After controlling the temperature setpoint, the amount of load shedding is given by Equation (20), where P DR k denotes the decrease of the HVAC load compared with the baseline situation. To measure the occupants' thermal discomfort, a discomfort function is introduced in Equation (21). The form of the discomfort function can be flexibly specified depending on the requirement on the temperature deviation. Let r k denote the unit revenue for the load shedding, then the optimal temperature control problem is described by Equation (22), where N k is the number of time intervals, T lower is the lower temperature limitation and T upper is the upper temperature limitation.
When the ramp rate of the indoor temperature is not constrained, the temperature setpoint only affects the loads of the current and the next time intervals. The decision variables are decoupled, and the optimization problem can be simplified into a series of optimization problems with a single decision variable, as shown in Equation (23). The physical explanation is that the proposed HVAC load model ignores buildings' longterm thermal inertia, so the influence of the current temperature only lasts until the next time interval. This decreases the performance of the model, but if the time interval is relatively large, such as half an hour, the errors are ignorable.
for k ∈ f1; 2; …; N k g We assume that the penalty is a quadratic function of the temperature deviation from the reference value, as shown in Equation (24). The choice of this function can be explained by the fact that people's thermal discomfort is usually small when the temperature deviation is rather small but increases rapidly when it exceeds a certain acceptable range. The smoothness of the quadratic function also simplifies the process of optimization:

| CASE STUDY
In this section, EnergyPlus is used to evaluate the performance of the proposed methods and to generate simulation data to train the proposed model. EnergyPlus is an open-source building energy simulation engine funded by the U.S. Department of Energy that has been widely applied to analyzing the energy consumption of HVAC systems. The test building model is an EnergyPlus example file whose geometric structure is shown in Figure 5. It is a single-story building located in Chicago with five rooms and a roof insulation layer. The run period of the simulation is from 1 June to 31 August, for 92 days in total. The weather data of Chicago is provided on the EnergyPlus website [33]. The simulation data are collected every hour, while the time step of the EnergyPlus simulation is 15 min. The simulation parameters are shown in Table 1.
We assume that the building participates in the real-time market and produces negawatts by load shedding [34]. Therefore, the electricity price in the real-time market is taken as the price signal for the load shedding. Hourly real-time price data from PJM for 1 June to 31 August are used in the simulation.

| Performance of HVAC load model
First, data are generated by running EnergyPlus simulation with the stepwise temperature setting strategy mentioned in Section 4.2. The simulation data of the total 92 days are split into the training set, the former 72 days, and the test set, the latter 20 days. The backfitting algorithm is used to find the parameters of the HVAC load model, and the number of iterations is set to be 20. Algorithms are programmed in the MATLAB R2019b and evaluated on a Lenovo ThinkPad P14s laptop with an Intel Core i7-10500U 1.8-GHz CPU and 16 GB of memory. The total computational time for finding the parameters is 27 s (the average of 10 individual experiments). The mean absolute percentage errors of the proposed model on the training set and the test set are 2.93% and 4.56%, respectively. A comparison between the forecast value and true value in the test set is depicted in Figure 6. The HVAC system is available between 6:00 and 20:00, and its electric power is zero during the rest of the day. For better display, only working hours are included in Figure 6. The results demonstrate the fitting and generalization capabilities of the proposed HVAC load model.

| Performance of optimal control strategy
To decide the coefficient c of the discomfort function in (5.2.1), several simulations are conducted beforehand. Keep the indoor temperature setting constant during the 3 month simulation period, then vary the indoor temperature setting and run the simulation multiple times with the same environmental data. Conduct linear regression to obtain how much HVAC load decreases with a unit increase in -219 temperature setting, which is denoted as a 0 . Calculate the average electricity price during the simulation period � r and multiply it by a 0 to obtain the numerical value of c. Note that the unit of c is not the same as a 0 � r. For the test building, the results are a 0 ¼ 0:645kW=℃, � r ¼ 24:3193$=MWh and c ¼ 0:0157$=ð℃Þ 2 . This step only aims to find a reasonable value for the coefficient c, which will be used to perform optimization. In practice, the coefficient c (or even the discomfort function) can be adjusted by the building operator, and a smaller value represents a greater willingness to reduce electricity cost at the expense of increased thermal discomfort. After determining the discomfort function, a simulation is conducted with the optimal control strategy. In fact, with the thermal discomfort penalty taking the quadratic form in Equation (24), the solution of Equation (23) can be directly written as Equation (25). The computational time for finding the optimal temperature setpoints is only 2e-6 s on average, which is rather ignorable.
To evaluate the control strategy, the simulation with the fixed indoor temperature setting is taken as the baseline. Taking the last 5 days of the simulation period for example, Figure 7 shows the control results, where (a) depicts the real-time price, (b) depicts the indoor temperature and (c) depicts the load shedding compared with the baseline. During the entire simulation period, the maximum and the minimum of the indoor temperature are 25:0℃ and 22:6℃ respectively. The cumulative cost saving compared with the baseline is chosen as the index of a strategy's performance. Figure 8 depicts the cumulative cost saving of the stepwise strategy and the proposed control strategy, where the total cost saving is the sum of the electricity cost saving and the discomfort cost saving. Compared with the stepwise strategy, the optimal control strategy saves more electricity cost with a relatively small increase in the discomfort cost. As for the cumulative total cost of the whole simulation period, the F I G U R E 7 Results of the optimal control strategy. (a) real-time price, (b) indoor temperature, (c) load shedding compared with fixed setting F I G U R E 8 Performance comparison 220 -HE ET AL. optimal control strategy saves 44.74% more than the stepwise strategy does.

| Discussion
The proposed optimal control strategy can adjust the temperature setting flexibly based on the electricity prices of the current and the next time step as well as the building's requirement for indoor temperature regulation. The performance of the control policy depends on the accuracy of the load forecast model because the impact of adjusting the temperature setpoint on the load is required to solve the optimization. The non-parametric components in the load model do not appear in the optimization, but without considering their effects the coefficients of the linear components cannot be accurately estimated.
The multiperiod thermal inertia of the building is ignored in the load model, and the relatively high performance confirms that this assumption is reasonable. This indicates that the effects of a building's thermal inertia usually last for a short period, specifically about an hour. As a result, precooling only happens when a sharp increase is expected in the next time step, and storing thermal energy in a building for a longer time appears impractical.
Because of the capacity limitation of the HVAC system, the actual indoor temperature does not always reach the setpoint. The change in temperature is not very dramatic, so the deviation between actual indoor temperature and the setpoint is ignored. When this deviation is significant, the control strategy should be modified; this is not included in the scope of this paper.

| CONCLUSION
A partially linear model is proposed to capture the electrothermal characteristics of building HVAC loads. The partially linear model combines a linear and a non-parametric model and thus has both explainability and flexibility. The model is also compatible with existing studies based on the assumption of a simplified RC model. The case study demonstrates the high accuracy of the load model. Based on the model, an optimal control strategy is developed to minimize building HVAC operating costs by providing negawatts in the real-time electricity market while considering the penalty of occupant discomfort. The case study shows that the proposed strategy can achieve a relatively high cost saving.

ACKNOWLEDGEMENT
This work was supported by the National Key R&D Programme of China under Grant No. 2020YFB0905900.