Abstract

This paper presents a greedy optimization algorithm for sampling design to calibrate WDS hydraulic model. The proposed approach starts from the existing sensors and sequentially adds one new sensor at each optimization simulation step. In each step, the algorithm tries to minimize the calibration prediction uncertainty. The new sensor is installed in the location where the uncertainty is greatest but also sensitive to other nodes. The robustness of the proposed approach is tested under different spatial and temporal demand distribution. We found that both the number of sensors and the perturbation ratio affect the calibration accuracy as defined by the average nodal pressure deviation itself and its variability. The plot of the calibration accuracy versus the number of sensors can reasonably guide the trade-off between model calibration accuracy and number of sensors placed or the cost. This proposed approach is superior in calibration accuracy and modeling efficiency when compared to the standard genetic algorithm (SGA) and Monte Carlo Sampling algorithm (MCS).

1. Introduction

Water distribution system (WDS) models are widely used by planners, water utility personnel, consultants, and many others engaged in the analysis, design, operation, and maintenance of water distribution networks [13]. WDS models can estimate flows and pressures in water distribution systems and further help understand contaminant transport and simple chemical interactions [47]. Thus, the model simulation can compensate for the lack of direct sensor measurements in the supervisory control and data acquisition (SCADA) system that collects real-time data of pipe flow, nodal pressure, and pump status [8]. The application of these models, however, depends on whether they are well calibrated in model parameters (e.g., the nodal consumer demand and the pipe roughness). In general, the model parameters are calculated to match observed data in acceptable levels [911], which are commonly field measurements collected using a set of sensors. Thus, the sampling design method is critical and important to the success of WDS model calibration.

The overall accuracy of a calibrated-model is generally improved with the increase of the number of sensors. However, the installation and maintenance cost increase as the number of sensors increases. A good trade-off between the number of sensors and the overall accuracy of the calibrated-model is desirable. Unnecessary cost should be avoided to collect the redundant data whose information contribution is already provided within other existing sensors [2].

Many methods have been proposed to determine how many and where the sensors should be installed [12, 13]. Most of them are involved in the use of Jacobian matrix (sensitive matrix) to minimize the calibration uncertainties caused by monitoring errors [1421]. These methods analyze the propagation of monitoring data uncertainties in the calibration process to minimize the uncertainty of the model parameters(roughness, nodal demand) or the model output (nodal pressure, and pipe flow), which can be classified under three different categories, namely, D-optimality criteria, A-optimality criteria, and V-optimality criteria [11]. D-optimality criteria and A-optimality criteria try to minimize the model parameter uncertainties while V- optimality criteria minimize the model output uncertainties. Particularly the A-optimality criteria minimize the variance matrix for the model parameter, whereas D-optimality criteria maximize the determinant of the inverse of the parameter variance matrix.

Bush et al. [15] proposed three methods for sampling design, namely, the max-sum, max-min, and weighted-sum methods, to minimize the parameter uncertainties. The three methods rank a potential sensor location based on the sensitivity matrix. In the max-sum method, the measurements are ranked on the basis of the summed absolute normalized sensitivity coefficients. The max-min method is also a ranking method to estimate all model parameters at once. For the weighted-sum method, the sensors are ranked by their contribution to the estimation parameters. Piller et al. [16] attempt to minimize the influence of measured error in the estimated parameter by selecting the rows of the Jacobian pseudoinverse with minimum matrix norms. Farley et al. [22] transform the sensitivity matrix into a 0-1 matrix and then calculate the sum of the elements of the columns of the 0-1matrix, selecting the largest element of the column to arrange the pressure sensors. Pérez et al. [21] obtain the 0-1 matrix in a similar way that uses an optimization algorithm to select the optimal combination of sensor points. Kapelan et al. [18] take the model accuracy and the cost into account by which two objective functions are formulated to optimize the sampling design through a multiobjective genetic algorithm. Kapelan et al. [19] use both the single-objective and multiobjective genetic algorithm to solve the sampling design problem.

However, all these methods, which are based on the analysis of Jacobian matrix, require estimating the calibration parameter values prior. This is difficult to achieve as model parameter can be estimated accurately only after the model calibration [20]. Generally, the sampling design is established before the calibration of the WDS model. Yet to solve the sampling design problem, WDS modeling needs to be conducted. More specifically, the Jacobian matrix and the prediction value must be generated through the WDS modeling. Unfortunately, the model parameter (nodal consumer demand and pipe roughness) of the WDS model has to be calibrated using sensor measurements. Thus, the sampling design problem resembles a “chicken and egg” dilemma [2]. To address this issue, a common practice is to generate the Jacobian matrix and the model prediction value that is referred to here as initial-model. These initial-model values can be just estimations (e.g., nodal consumer demand and pipe roughness) based on existing data [9, 15, 18, 19]. However, the parameters of the initial-model may be far away from the real water distribution system or the ultimate global optimization state. Consequently, model prediction of the initial-model can differ significantly from the true values, for which the sampling design methods based on the initial-model may be unreliable. Therefore, robustness of the optimization method is at the center of technical consideration. To examine the model fidelity, performance of a sampling design method is normally tested under different initial-models with varied spatial and temporal demand distributions.

This paper proposes a greedy algorithm to optimize the sampling design for the calibration of WDS hydraulic model. Here we are focused on the calibration of consumer demands; the proposed approach can be used for calibration of other parameters (e.g., pipe roughness). Generally, the greedy algorithm starts from the existing sensors and sequentially adds one new sensor at each step. The new sensor is selected from the candidate sensor locations that have the ability to improve the overall calibration accuracy. The robustness of the proposed approach is tested by 10000 runs of varied demand distributions. Its performance (the calibration accuracy and efficiency) is compared with two global optimal algorithms, namely, the standard genetic algorithm (SGA) and Monte Carlo Sampling (MC).

2. Methodology

An outline of the content and logical structure of this paper is given in Figure 1, with details presented in the subsequent sections.

2.1. Definition of WDS Calibration Problem

WDS calibration is a process of comparing model prediction value (nodal pressure and pipe flow) with measured value and making the appropriate adjustments of the model parameter so that both results agree. The water demand calibration objective is typically to adjust the water demands within their feasible domains in the direction of reducing the calibration objective [23]. Typically, the water demand calibration aims to minimize the square difference between the measured value and the prediction value [24] in order to find the consumer demand of the calibrated hydraulic model, which is referred to as calibrated-model in this paper. The objective function can be stated as

subject towhere is the objective function; is the vector of consumer demand, which is the decision variable to be solved; and are the measured and prediction value at node (pipe) ; is the number of measurements; are the weighting factors for the measurements; is the system of nonlinear equations describing the hydraulic steady state of flows and pressure in a water distribution network, which includes mass continuity and energy conservation equations.

The calibration problem of the WDS is an inverse problem [25, 26]. We rely on monitoring data to estimate real model parameters. In practice, the number of the sensors is less than the number of nodal consumer demand. Thus, model calibration is an ill-conditioning problem [2729], which will result IN the parameter uncertainty [30]. For the nodes (pipes) with the sensors, we can minimize (1) so that the prediction value of the calibrated-model fits the measured value, while for nodes with no sensors, the prediction value may deviate from true values [29].

In this paper, (1) is minimized by the calibration algorithm proposed by Cheng et al. [24]. This algorithm minimizes (1) by singular value decomposition to calculate search director, which is more reliable and effective than the finite-difference method [24].

2.2. Sampling Design Optimization Approach

The sampling design optimization approach is to find the optimal set of network locations for the sensors to collect data for the calibration of WDS model [19, 20]. As mentioned previously, the sampling design approaches are typically based on the analysis of the initial-model with the calibration parameter values estimated prior. It treats N number of prediction values (nodal pressure and pipe flow) selected from the initial-model as the measured value [31], referred to as pseudo-measured value. The pseudo-measured value is used to calibrate the WDS model by minimizing (1) to generate a calibrated-model. Typically, the calibrated-model is deviated from the initial-model as the pseudo-measurements provide only partial information about the initial-model. The closer the calibrated-model is to the initial-model, the more effective the information provided by the pseudo-measurements is. In this paper, the calibration prediction uncertainty is evaluated by the predictions deviation between calibrated-model and initial-model, which are stated in (3). Thus, the objective function of the sampling design is to minimize the average deviation, which is stated in (4).where is the vectors of nodal pressure deviation, which is referred to as calibration deviation vector; is the pressure prediction value of node in the initial-model; is the pressure prediction value of node in the calibrated-model; is the number of nodes.

Equation (4) will be minimized when the number of sensors is increased, since the nodal pressure deviation is reduced with the increase of the number of sensors. On the one hand, the prediction deviation corresponding to the sensors can be minimized by (1). On the other hand, the pressure deviations of the nodes that are sensitive to the measurements are also reduced. Thus, the calibration accuracy improvement consists of two aspects: (1) the prediction deviations where the sensor located can be reduced and (2) the prediction deviations where there are no sensors installed but they are sensitive to the sensor that can also be reduced. The sampling design optimization needs to maximize these two improvements, which are quantified as follows.

2.2.1. The Calibration Accuracy Improvement at the Pressure Sensor Node

Assuming that some sensors already exist in the network, we can generate the calibrated-model by minimizing (1). Then is calculated by (3). In order to improve the calibration accuracy, we need to increase the number of sensors. To minimize (4), a new sensor should be installed at the node (pipe) which has the maximum element value . This is because we can minimize the calibration deviation by installing the sensor at the location with the greatest deviation. For example, assuming that node will be selected to install a new pressure sensor, the deviation of node will reduce to zero by minimizing (1).

This approach can effectively reduce the deviation of the sensor nodes. However, the number of sensors is far less than the number of nodes without sensors. The pressure deviation for these nodes without sensors should be also reduced.

2.2.2. The Calibration Accuracy Improvement of All the Nodes

The pressure deviation of the nodes which are sensitive to the sensors can be also reduced. Thus, sensors are installed in places which have the maximum sensitivity on pressure variation on it.. The pressure sensitivity matrix () reflects the sensitivity relationship between the pressures of each node. means the effect of the pressure variation at node on the pressure variation at node . For example, assuming that the pressure variation at node is 1, the pressure variation at node will be . Farley et al. [22] proposed a method to get the pressure sensitive matrix by giving an extra flow in turn on each node in the initial-model as follows:(1)Run EPANet to simulate the initial-model and obtain nodal pressure .(2)Give an extra flow at node and obtain nodal pressure .(3)Calculate the elements of : , being the absolute value of variation of pressure at node after adding an extra flow at node . is the pressure of node in initial-model and is the pressure of node after adding extra flow at node .(4)Reset nodal demand to original value and add another extra flow to next node.

This work generates the pressure sensitive matrix :

As mentioned previously, describes the effect of the pressure variation at node on the pressure variation at node . The sensors with different locations are not equally sensitive to the WDS hydraulic system [21, 32, 33] and the overall pressure sensitivity of the node can be quantified by . Assuming that a new pressure sensor will be installed at node , the calibration accuracy improvement at node itself will be . Correspondingly, the calibration deviation of node will be reduced by . Therefore the pressure deviation of all the nodes are expected to be reduced by . The pressure sensor should be then installed at the node with the maximum of , where is the element of which is referred to as “pressure weight coefficient vector”.

2.3. The Process of the Sampling Design

The proposed sampling design approach starts from the existing sensors and selects one new sensor at each optimization simulation step. In most cases, the WDS networks are already equipped with a certain number of monitoring equipment, especially for the water station or the pump station. In the absence of sensors, the first sensor should be started at the water stations (reservoir). The pseudo-measured value and the pressure weight coefficient vector () are estimated a priori based on the simulation the initial-model. In each optimization simulation step, the pseudo-measured value corresponding to the existing sensors are used to generate the calibrated-model by solving (1). In this paper (1) is minimized by the calibration algorithm proposed by Cheng et al. [24]. Then the calibration deviation vector is obtained by calculating the deviation between initial-model and calibrated-model. Thus, a new sensor is selected at the node with the maximum of . For the next step, the new added sensor combined with other sensors is treated as the existing sensors to calibrate the WDS modal, calculate the deviation, and then select next new sensor. Figure 2 shows the process of the sample design for the pressure sensor.

The sampling design for the flow sensors is carried out after completing the pressure sampling design. This is because the nodal pressure is much sensitive than the pipe flow rate for water demand calibration. The pressure measurements can improve the pressure calibration accuracy of other nodes and also have a considerable impact on the calibration of pipe flow rates. Therefore, priority placement of pressure sensors is conducive to the overall improvement of the calibration accuracy.

The flow sensors are also selected one by one. In each step, the pipe with biggest flow deviation () will be selected to place the new flow sensor. where are the vectors of pipe flow deviation; is the flow prediction value of pipe in the initial-model; is the flow prediction value of pipe in the calibrated-model; is the number of pipes.

The greedy algorithm relies on initial-model to derive the pseudo-measured value and the pressure weight coefficient vector (). Typically, initial-model needs to be as accurate as possible. The nodal consumer demand of the initial-model can be estimated a priori by experience or billings.

3. Result and Discussion

3.1. Application to a Real Network

The proposed approach is applied to a real network in Figure 3. This network contains 3 water supply stations, 491 nodes, and 640 pipes. Water outflows of the water supply stations are monitored, as the existing flow sensors. The total consumer demand is the sum of the outflows of the water stations. The consumer demands of initial-model are assigned by allocating the total consumer demand to the nodes which are weighted by the billing records. Figure 3 also shows the final optimized placement of 20 pressure sensors and 10 flow sensors with step-by-step increment.

As mentioned in the Methodology, only one new sensor is added at each step. The new sensor should not only improve the calibration accuracy of itself but also improve the calibration accuracy of other nodes. This can be illustrated by Figure 4. Figure 4(a) shows the pressure deviation of all the nodes when the number of sensors is 10. One node (the black tetragonal node) will be selected as the new sensor in the next step. The circle nodes are found to be sensitive to this new sensor node according to the pressure sensitive matrix. Figure 4(b) shows the pressure deviation of the nodes when the new sensor has been installed. Comparing Figure 4(b) with Figure 4(a), we can see that pressure deviations of the sensor node (the black tetragonal node) decrease from 5.29 m to 0 m. Besides, the pressure deviation of the circle nodes also decreases obviously, but less than that of the sensor node.

Figure 5 shows the average pressure deviation versus the number of sensors. The average pressure deviation is calculated by (4). Generally, the deviation decreases as the number of sensors increases. As shown in Figure 5, the average pressure deviation is 1.16 m for 2 pressure sensors and 0 flow sensor. When the number of pressure sensors reaches up to 30, the corresponding average pressure deviation becomes stabilized at 0.19 m. The average pressure deviation also decreases as the number of flow sensors increases. As the number of flow sensors increases from 0 to 15 consistently with 2 pressure sensors, the average pressure deviation decreases from 1.16 m to 0.27 m.

In general, the more the number of sensors there is, the more likely it is to obtain good calibration accuracy. However, we can also find that the pressure deviations do not always decrease with the increase of the number of sensors when there are only a few sensors. For example, one new pressure sensor is added to the existing 3 pressure sensors without flow sensor; the average pressure deviation slightly increases from 1.1 m to 1.17 m. A similar phenomenon also happens in 2-3 pressure sensors with 5 flow sensors. A possible explanation could be that fewer sensors lead to severe ill problem for the calibration problem. Many possible solutions could be accepted, and the calibration method [24] chooses one solution from many possible solutions according to their algorithm bias. Thus, the improvement of the calibration accuracy by adding a new sensor is not ensured, especially when the number of sensors is small.

It is interesting to note that two different performance-changing stages can be observed from Figure 5. For the first stage with the number of pressure sensors within 16, it shows a rapid decrease of the average pressure deviation; then followed by the second stage with the number of pressure sensors exceeding 16 the average pressure deviation decreases in a moderate rate. The effectiveness-cost ratio of sensors placement diminishes with increasing sensors. This may be attributed to the sensor selection strategy. In the first stage, the number of improved nodes increases as the number of sensors increases. In the second phase, the nodes which are sensitive to the new sensor may be also sensitive to the existing sensors. Thus, the contribution of the new sensor is already contained within the existing sensors.

3.2. The Robustness Tests

The proposed approach is based on the initial-model and the calibration accuracy is illustrated by analyzing the deviation between calibrated-model and initial-model. A fundamental issue regarding calibration accuracy is that the true demand distribution is often unknown and the sensitive matrix cannot fully and accurately reflect the sensitivity relationship between nodal pressures, but can be an approximation. Moreover, the sensors placement scheme should be applied to various demand distributions, which may be much different from the initial-model. We should evaluate the performance of the sensor placement scheme generated by the proposed approach under other demand distributions. Thus, robustness tests are designed to address this issue. 10000 varied demand distributions are generated by giving a random perturbation to the consumer demand of initial-model. The demand distributions are generated according to (9).where is the consumer demands of node in the initial-model and is the consumer demand of node after giving a random perturbation to . represents the uniform distribution. is the perturbation ratio; is a random value which obeys the uniform distribution with .

The perturbation ratio () represents the deviation degree of the demand distribution from initial-model. For example, for =50%, it means that the consumer demand of the demand distribution has 50% deviation to initial-model. In this paper, different values were tested. However, is limited to within 100%, because a larger value can cause unreasonable demand distributions, such as a small pipe with a large flow and a large pipe with a small flow.

The sensors placement scheme generated from the initial-model is used to calibrate the model with abovementioned demand distributions. The calibration accuracies for those models are analyzed. The process of the robustness test is shown in Figure 6. The average pressure deviation for 10000 different demand distributions with corresponding perturbation ratios is illustrated in Figure 7. The proposed approach shows a good robustness since the average pressure deviation under different demand distributions is consistent with the deviation of initial-model (=0%) as the number of pressure sensors increases. Two different performance-changing stages separated by 16 pressure sensors can be observed from the different demand distributions, similar to the case before. The average pressure deviation decreases rapidly in the first stage and decreases in a moderate rate in the second stage. This means that we can use the initial-model to determine the appropriate number of sensors, and it will not be significantly affected by the actual water demand distribution.

It is worth mentioning that the calibration accuracy of the sensor scheme is slightly worse with the increase of perturbation ratio (). For example, when =10%, the average pressure deviation is almost the same as that of initial-model. The average pressure deviation increases about 0.12 m compared to initial-model when =100%. This indicates that the sensor placement scheme obtained from the initial-model has the ability to calibrate the hydraulic model under other demand distributions effectively.

The variability of the average pressure deviation due to different demand distributions can be quantified by the interval between the upper and the lower bounds of the 95% confidence interval. Figure 7 shows that the variability increases with the increase of the perturbation ratio. This is consistent with the average pressure deviation which also increases with the increase of perturbation ratio. One possible reason is that the number of nodes which are sensitive to the sensors is determined by the sensitivity matrix. But with the increase of , the sensitivity matrix of the initial-model may be no longer applicable. For those nodes that the accuracy cannot be improved, their consumer demand fluctuates wildly with the increase of , which in turn leads to large variability in the pressure deviation of these nodes. This validates the previous conclusion that deviation degree of the initial-model from the true demand distribution has slightly influence on the calibration accuracy.

The average pressure deviation variability versus sensors number curves also has the characteristic of two different stages which are separated by 16 pressure sensors. The deviation variability decreases rapidly in the first stage, and it almost no longer decreases in the second phase. For example, when =100%, the variability interval of 95% confidence interval decreases from 1.11 m to 0.37 m as the number of pressure sensors increases from 1 to 16, while it decreases from 0.3 m to 0.25 m when the number of pressure sensors increases from 16 to 30. All of this validates a conclusion that both the deviation and the variability can be reduced by increasing the number of sensors in the first stage. As for the second stage, it is not economical to increase the number of sensors since installing of new sensors has a bad cost-efficiency ratio.

3.3. The Accuracy Compared with Global Optimal Algorithms

As mentioned previously, the proposed greedy approach cannot guarantee the global optimum, since the sensors are selected one by one. In order to test the reliability of this approach, the calibration accuracy of the proposed approach is compared with two different global optimal algorithms, namely, the standard genetic algorithm and (SGA) Monte Carlo Sampling algorithm (MCS).

The SGA was introduced by Holland [34]. The SGA relies on three fundamental operators, namely, selection, crossover, and mutation. These operators are mainly driven by three parameters: the crossover rate (), mutation rate (), and the population size (). In this paper, and are obtained by trial-and-error analysis. For this case, = 0.8 and =0.15. An initial population with 100 initial solutions is generated randomly, corresponding to = 100. Each solution consists of elements, corresponding to sensors and each element is the index of the selected sensor node. The maximum number of iterations is 1000. The objective function is to minimize (4), which is the same as the proposed approach.

The MCS generates the sensors’ location by random sampling. In this study, 10000 sensor placement schemes are generated randomly, and only the best scheme which has the minimum deviation is selected.

A comparison of the performance between the proposed approach and the two global optimization algorithms (SGA and MCS) is shown in Figure 8. The overall pressure deviation for the three approaches decreases with the increase of sensor numbers. The proposed approach performs better when the number of sensors is large. As can be seen from Figure 8, when the number of sensor is within 9, the proposed approach performs worse than other two algorithms. When the number of sensors ranges from 9 to 15, the proposed approach is still slightly worse than SGA but much better than MCS. As the number of sensors exceeds 15, the proposed approach performs better than both SGA and MCS. This implies that the proposed approach is more sufficient for a large number of sensors' placements. One possible explanation is that the search spaces of the SGA and MCS increase with the increase of the number of sensors. When the number of sensors is , the number of candidate solutions is for SGA and MCS and is the number of the nodes of WDS model. For the greedy algorithm, the number of candidate solutions is since only one sensor is selected from the nodes every step. For small number of sensors, SGA and MCS can find the global optimal which is similar to go through all the potential solutions, while with the increase of the number of sensors, it becomes more and more difficult for SGA and MCS to find the optimal solution. Thus, the proposed approach performs better than other two algorithms for the large number of sensors. However, even the global optimization approaches still do not guarantee that pressures sensors perform better than pressure sensors. In this aspect, the proposed approach is better than these global optimization approaches when the number of sensors is large.

The proposed approach shows a superior performance in time-consuming over the other two approaches. Figure 9 shows that the time-consuming of the three proposed approaches is 101, 107, and 108 ms for 1 pressure sensor. When the number of sensors reaches up to 30, the time-consuming of the three approaches is 105, 107, and 109 ms.

4. Conclusions

This paper proposed a greedy sampling design algorithm to select the sensor location with the greatest improvement for optimal calibration accuracy. This improvement is quantified by the calibration deviation vector and the pressure weight coefficient vector. Only one pressure sensor is selected at each algorithm step. This sampling design algorithm can be based on any model calibration method. Its final result may be different relied on which model calibration method is chosen. The sampling design algorithm should be used in conjunction with the corresponding model calibration methods.

The use of greedy sampling design approach may solve the “chicken and egg” dilemma commonly encountered in the sensor placement for model calibration. Through a real-world case study, it is found that, by using the greedy sampling design algorithm, the ability to calibrate WDS model with wildly changed demand distributions has improved with the sensors number no more than 20. The calibration accuracy will slightly decrease as the deviation between initial-model and true model increases.

In comparison with those of global optimization algorithms (i.e., SGA and MCS), the proposed approach performs better than SGA and MCS when the number of sensors is large, while it performs slightly worse than SGA and MCS when the number of sensors is small. This is because the search space of the proposed greedy algorithm decreases for every algorithm step, while the search spaces of SGA and MCS expand exponentially as the number of sensors increased.

The relationship between the number of sensors and the calibration accuracy is also quantified. The proposed approach is applied to a real-world network, and two different performance-changing stages can be observed. The average pressure deviation decreases rapidly in the first stages and decreases in a moderate way in the second stage. Moreover, the average pressure deviation variability versus sensors number curves also has the characteristic of two different stages. This can reasonably guide the trade-off between model calibration accuracy and number of sensors placed or the cost.

Notations

The following symbols are used in this paper:
:Average pressure deviation
:Crossover rate
:Nonlinear equations of the WDS network
:Vector of nodal pressure
:Absolute value of variation of nodal pressure
:Prediction value at node (pipe)
:The calibration objective function
:Mutation rate
:Number of nodes
:The number of pipes
:Number of samples
:Population size
:Vector of consumer demand
:Random value
:Perturbation ratio
:Vector of pipe flow
:Pressure sensitive matrix
:The uniform distribution
:Pressure weight coefficient vector
:Weighting factors for the measurements
:Measured value at node (pipe)
: Calibration deviation vector for nodal pressure
: Calibration deviation vector for pipe flow.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure

Although not directly funded by US EPA, this paper has also been subjected to the Agency’s administrative review and has been approved for external publication. Any opinions expressed in this paper are those of the authors and do not necessarily reflect the views of the Agency; therefore, no official endorsement should be inferred.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The present research is funded by the National Key Research and Development Program of China (No. 2016YFC0400600), Science and Technology Program of Zhejiang Province (Nos. 2017C33174 and 2015C33007), and the National Natural Science Foundation of China (Nos. 51208457 and 51478417).