Optimal Combinations of Non-Sequential Regressors for ARX-Based Typhoon Inundation Forecast Models Considering Multiple Objectives

Inundation forecast models with non-sequential regressors are advantageous in efficiency due to their rather fewer input variables required to be processed. This type of model is nevertheless rare mainly because of the difficulty in finding the proper combination of regressors for the model to perform accurate prediction. A novel methodology is proposed in this study to tackle the problem. The approach involves integrating a Multi-Objective Genetic Algorithm (MOGA) with forecast models based on ARX (Auto-Regressive model with eXogenous inputs) to transfer the search for the optimal combination of non-sequential regressors into an optimization problem. An innovative approach to codifying any combinations of model regressors into binary strings is developed and employed in MOGA. The Pareto optimal sets of three types of models including linear ARX (LARX), nonlinear ARX with Wavelet function (NLARX-W), and nonlinear ARX with Sigmoid function (NLARX-S) are searched for by the proposed methodology. The results show that the optimal models acquired through this approach have good inundation forecasting capabilities in every aspect in terms of accuracy, time shift error, and error distribution.


Introduction
Typhoon is a common weather phenomenon in subtropical areas and usually occurs between July and October of each year.Heavy rainfall brought by the typhoon during the event usually results in serious inundation problems in low-lying areas, which not only causes property loss to the local population but also threatens their safety.Due to restrictions on engineering funding, structural protective measures are constrained by the designed limits.Once the typhoon scale exceeds the designed protective limit, people must rely on nonstructural means for disaster relief during the event, such as evacuating people from areas in potentially high flooding risk.Among nonstructural measures the accurate forecast of the inundation level in the areas within the next several hours is a critical factor in the decision-making and planning of disaster relief actions.
Relevant studies on inundation forecast technology are quite ample and can generally be divided into either the numerical simulation or the black-box modeling.The numerical simulation is based on theoretical deduction through the understanding of the mechanism from rainfall to inundation.The advantage of this type of method is the completed support from the theoretical basis for the physical mechanism of inundation.The simulation result often has a high degree of accuracy, which renders the method a powerful tool of inundation forecast.However, the disadvantage of this method is its high demand for computing resources and CPU time, which makes it unsuitable to provide the real-time forecast required in the quick disaster prevention and rescue response during the typhoon attack.On the other hand, the black-box modeling relies on a different approach by deeming the process from rainfall to inundation as a black box.It does not delve into the internal physical mechanism but instead analyzes the input and output data of the system to simulate the relationship between them.These types of models cannot explain the physical mechanism involved in the system, but they can correctly and effectively simulate the response of the system, and the computing speed is faster than numerical models [1].These practical benefits render black-box modeling a suitable forecasting tool for decision making and rescue planning during the typhoon period.
Abundant studies with regard to black-box modeling for inundation forecasting can be found in literature, for example, Liong et al. [2] developed a river stage forecasting model based on an artificial neural network (ANN) and yielded a very high degree of prediction accuracy even for up to seven lead days.Campolo et al. [3] developed a flood forecasting model based on ANN that exploits real-time information available for the basin of the River Arno to predict the basin's water level evolution.Keskin et al. [4] proposed a flow prediction method based on an adaptive neural-based fuzzy inference system (ANFIS) coupled with stochastic hydrological models.Shu and Ouarda [5] proposed a methodology using ANFIS for flood quantile estimation at ungauged sites and demonstrated that the ANFIS approach has a much better generalization capability than other alternatives.Kia et al. [6] develop a flood model based on ANN using various flood causative factors in conjunction with geographic information system (GIS) to model flood-prone areas in southern Malaysia.Lin et al. [7] proposed a real-time regional forecasting model to yield 1-to 3-h lead time inundation maps based on K-means cluster analysis incorporated with support vector machine (SVM).Tehrany et al. [8] proposed a methodology for flood susceptibility mapping by combining SVM and weights-of-evidence (WoE) models and demonstrated that the ensemble method outperforms the individual methods.Del Giudice et al. [9] developed a methodology that formulated models with increasing detail and flexibility, describing their systematic deviations using an autoregressive bias process.Chang and Tsai [10] proposed a spatial-temporal lumping of radar rainfall for modeling inflow forecasts to mitigate time-lag problems and improve forecasting accuracy.
In view of the above literature, most of the approaches employ sequential data as model inputs.Less has been explored for forecasting models with non-sequential data inputs.This type of model is efficient because there are rather fewer inputs required to be processed.The challenge for these models, however, lies in selecting the appropriate combination of non-sequential variables to be used in the inputs.This study aims to propose a methodology to tackle this difficulty.The approach proceeds by integrating a Multi-Objective Genetic Algorithm (MOGA) with models based on ARX (Auto-Regressive models with eXogenous inputs) to search for the optimal combination of non-sequential regressors for model inputs.Three types of ARX-based models are tested by the proposed methodology, including linear ARX (LARX), nonlinear ARX with Wavelet function (NLARX-W), and nonlinear ARX with Sigmoid function (NLARX-S).The models are assessed by a number of indexes to examine their performance on various aspects, and the characteristics of the models selected from the Pareto optimal sets located by MOGA with the best performance in each index are compared and discussed.The remainder of this paper is arranged as follows: Section 2 illustrates the environmental background of the study area, the ARX models, and the optimal models obtained by MOGA; Section 3 discusses the characteristics of the acquired optimal models, and finally the conclusions are drawn in Section 4 based on the findings.

Study Area
Yilan County (Figure 1), located in northeastern Taiwan, is known for its rainy weather.With about 200 rainy days per year, the annual average precipitation is around 2000-2500 mm.The terrain is surrounded by mountains in the west and coastline in the east.It is frequently hit by typhoons in the summer and autumn of each year.On average, two to three typhoons hit Taiwan every year and among which 45% land in Yilan County [11].The heavy rainfall during the typhoon attacking Water 2017, 9, 519 3 of 15 period often causes severe inundation in the low-lying areas.Among these inundation-prone areas, the situation in Donsan area (Figure 1) is most extreme.The drainage area of Donsan is about 34 km 2 .Ground level in the area is very low and flat.The ground elevation in more than one-third of the area is below 2 m, as seen in Figure 1.During typhoon invasion, these low grounds are often flooded, causing heavy damages and tremendous losses to the properties.
Water 2017, 9, 519 3 of 15 situation in Donsan area (Figure 1) is most extreme.The drainage area of Donsan is about 34 km 2 .Ground level in the area is very low and flat.The ground elevation in more than one-third of the area is below 2 m, as seen in Figure 1.During typhoon invasion, these low grounds are often flooded, causing heavy damages and tremendous losses to the properties.The high tendency in severe flooding within the area has raised an urgent need for effective disaster preventive measures.In order to monitor the local inundation condition during the typhoon period, Taiwan Water Resources Agency established a surveillance network for the area in 2011.It includes three water-level gauging stations, as well as a data transmission system that can receive the precipitation observations of the Quantitative Precipitation Estimation and Segregation Using Multiple Sensor (QPESUMS, [12]).QPESUMS consists of eight Doppler radar stations that each scans an area with a radius of approximately 230 km.The system provides quantitative precipitation estimation by integrating observations from weather radars and rainfall readings from 406 automatic gauges and 45 ground stations in Taiwan.QPESUMS also provides rainfall forecasts by tracking and extrapolating the movement paths of storm cells based on radar readings.During the typhoon period, QPESUMS delivers the rainfall data in a frequency of every 10 min through a network connection to the surveillance system.The water level gauging stations return the local inundation water-level data with the same frequency through radio transmission.
Since the monitoring system of Donsan area was established, it has collected data of 10 typhoons over the years.Figure 2 shows the recorded water level of each station and QPESUMS rainfall data during Typhoon Trami in 2013.Detailed information of each typhoon event is listed in Table 1.These data not only provide the local rainfall and water level information during the typhoon period but can also be used for the establishment of water-level forecast models.In the present study, the WG2 gauging site is selected as the study object to establish the water level forecast model.The high tendency in severe flooding within the area has raised an urgent need for effective disaster preventive measures.In order to monitor the local inundation condition during the typhoon period, Taiwan Water Resources Agency established a surveillance network for the area in 2011.It includes three water-level gauging stations, as well as a data transmission system that can receive the precipitation observations of the Quantitative Precipitation Estimation and Segregation Using Multiple Sensor (QPESUMS, [12]).QPESUMS consists of eight Doppler radar stations that each scans an area with a radius of approximately 230 km.The system provides quantitative precipitation estimation by integrating observations from weather radars and rainfall readings from 406 automatic gauges and 45 ground stations in Taiwan.QPESUMS also provides rainfall forecasts by tracking and extrapolating the movement paths of storm cells based on radar readings.During the typhoon period, QPESUMS delivers the rainfall data in a frequency of every 10 min through a network connection to the surveillance system.The water level gauging stations return the local inundation water-level data with the same frequency through radio transmission.
Since the monitoring system of Donsan area was established, it has collected data of 10 typhoons over the years.Figure 2 shows the recorded water level of each station and QPESUMS rainfall data during Typhoon Trami in 2013.Detailed information of each typhoon event is listed in Table 1.These data not only provide the local rainfall and water level information during the typhoon period but can also be used for the establishment of water-level forecast models.In the present study, the WG2 gauging site is selected as the study object to establish the water level forecast model.situation in Donsan area (Figure 1) is most extreme.The drainage area of Donsan is about 34 km 2 .Ground level in the area is very low and flat.The ground elevation in more than one-third of the area is below 2 m, as seen in Figure 1.During typhoon invasion, these low grounds are often flooded, causing heavy damages and tremendous losses to the properties.The high tendency in severe flooding within the area has raised an urgent need for effective disaster preventive measures.In order to monitor the local inundation condition during the typhoon period, Taiwan Water Resources Agency established a surveillance network for the area in 2011.It includes three water-level gauging stations, as well as a data transmission system that can receive the precipitation observations of the Quantitative Precipitation Estimation and Segregation Using Multiple Sensor (QPESUMS, [12]).QPESUMS consists of eight Doppler radar stations that each scans an area with a radius of approximately 230 km.The system provides quantitative precipitation estimation by integrating observations from weather radars and rainfall readings from 406 automatic gauges and 45 ground stations in Taiwan.QPESUMS also provides rainfall forecasts by tracking and extrapolating the movement paths of storm cells based on radar readings.During the typhoon period, QPESUMS delivers the rainfall data in a frequency of every 10 min through a network connection to the surveillance system.The water level gauging stations return the local inundation water-level data with the same frequency through radio transmission.
Since the monitoring system of Donsan area was established, it has collected data of 10 typhoons over the years.Figure 2 shows the recorded water level of each station and QPESUMS rainfall data during Typhoon Trami in 2013.Detailed information of each typhoon event is listed in Table 1.These data not only provide the local rainfall and water level information during the typhoon period but can also be used for the establishment of water-level forecast models.In the present study, the WG2 gauging site is selected as the study object to establish the water level forecast model.

Model Construction
The present study adopts ARX to build the water-level forecast model for the quick calculation feature.ARX is a kind of black-box model which can be divided into the linear model and nonlinear model according to the relationship between the input and output.The nonlinear model can be further subdivided into various types according to the different nonlinear functions applied.The present study respectively aims at linear ARX, as well as ARX with Wavelet function and ARX with Sigmoid function in the nonlinear ARX, to discuss the inundation forecast performance of these three types of models.
With regard to emergency response to the typhoon, the water level is a greater concern than runoff.Therefore, the purpose of the forecast model is to establish the relationship between rainfall and water level, rather than the more commonly seen rainfall and runoff relationship.Although few studies have focused on the relationship between rainfall and water level, relevant works can be found in the literature [13][14][15].

Linear ARX
Linear ARX (LARX) is an extension of the AR model [16] in time series analysis by adding the influence of other exogenous inputs.The equation is as follows: in which H represents the system output; R is the exogenous input; n a and n b are the number of terms of H and R, respectively; n k is the time lag of R; e is white noise; a 1 through a n a and b 1 through b n b are the coefficient of each term, respectively.
are called the regressors which take in the known data to forecast H(t + 1).In LARX, the relationship between H(t) and the regressors is linear, as shown in Equation (1).In this study, H represents the water level at the site of WG2 in Donsan area, while R is the local rainfall data provided by QPESUMS.

Nonlinear ARX
Nonlinear ARX extends linear ARX by using a nonlinear relationship between the forecast and the regressors, as shown in the equation below: where ψ is the nonlinearity estimator which has the following form: in which x is the row vector consists of the regressors, n is the number of terms of nonlinearity estimator, α k and β k are the coefficients of each term, and γ k is the mean value of the regressor vector.κ is the nonlinear unit, which can adopt different nonlinear functions.
Two different types of nonlinear ARX models are employed in the present study, viz. the nonlinear ARX with Wavelet function (NLARX-W) and the nonlinear ARX with Sigmoid function (NLARX-S).In NLARX-W, the nonlinear unit κ adopts the Wavelet function [17], which has the following formula: where dim represents the dimension of the vector s; e is the exponential function.
On the other hand, in NLARX-S, κ adopts the Sigmoid function as shown in the equation below: The model structure of LARX, NLARX-W, and NLARX-S is determined by the selection of regressors which are to be determined later on through the optimization process.

Rainfall Data Analysis
To analyze the relationship between water level and rainfall, correlation analysis on typhoon data records was carried out.The definition of the correlation coefficient (CC) is shown below [18]: in which cov is the covariance between variables x and y; σ x and σ y are the standard deviations of x and y, respectively; n is the quantity of the data points.CC ranges between −1 and 1, where −1 represents the perfect negative correlation between x and y, +1 means the perfect positive correlation, and 0 means that there is no correlation between x and y.
The rainfall data provided by QPESUMS are 10-min rainfall increments.However, it has been reported by Ouyang [19] that the variation in water level is slower than the variation in rainfall, and often a certain amount of rainfall accumulation is required for the water to build up.To search for the most correlated cumulative rainfall to the water level at the study site, the procedure of Ouyang [19] was adopted here.The moving cumulative rainfall data of the typhoon events with accumulation duration ranging from 1 h to 30 h were first constructed.The correlations of each cumulative rainfall data set and the water level data were then calculated.The results are as shown in Figure 3.The circle dots in the figure represent the average CC of the typhoon events, and the error bars denote the CC distribution of the events.As shown in the figure, the average CC is first gradually increasing along with an increase in the accumulation duration, and then followed by a gradual recession with a further increase in the accumulation duration.As shown in the figure, the maximum average CC appears when the accumulation duration is about 18 h.This indicates that the water level at the study site is most correlated with the 18 h cumulative rainfall which shall be adopted as the input data of the ARX models.
where dim represents the dimension of the vector ; is the exponential function.
On the other hand, in NLARX-S, adopts the Sigmoid function as shown in the equation below: The model structure of LARX, NLARX-W, and NLARX-S is determined by the selection of regressors which are to be determined later on through the optimization process.

Rainfall Data Analysis
To analyze the relationship between water level and rainfall, correlation analysis on typhoon data records was carried out.The definition of the correlation coefficient (CC) is shown below [18]: in which cov is the covariance between variables and ; and are the standard deviations of and , respectively; is the quantity of the data points.CC ranges between −1 and 1, where −1 represents the perfect negative correlation between and , +1 means the perfect positive correlation, and 0 means that there is no correlation between and .
The rainfall data provided by QPESUMS are 10-min rainfall increments.However, it has been reported by Ouyang [19] that the variation in water level is slower than the variation in rainfall, and often a certain amount of rainfall accumulation is required for the water to build up.To search for the most correlated cumulative rainfall to the water level at the study site, the procedure of Ouyang [19] was adopted here.The moving cumulative rainfall data of the typhoon events with accumulation duration ranging from 1 h to 30 h were first constructed.The correlations of each cumulative rainfall data set and the water level data were then calculated.The results are as shown in Figure 3.The circle dots in the figure represent the average CC of the typhoon events, and the error bars denote the CC distribution of the events.As shown in the figure, the average CC is first gradually increasing along with an increase in the accumulation duration, and then followed by a gradual recession with a further increase in the accumulation duration.As shown in the figure, the maximum average CC appears when the accumulation duration is about 18 h.This indicates that the water level at the study site is most correlated with the 18 h cumulative rainfall which shall be adopted as the input data of the ARX models.

Regressors
The ARX simulation depends on the selection of regressors.With different combinations of regressors, the performance of the model varies.In the selection of model regressors, two methods are commonly utilized in the literature.The first one is Sequential Time Series (see for example, [20][21][22][23]) where a sequence of time series data from the current time to a certain time before is used as model regressors, for example, R(t), R(t − 1), . . ., R(t − m), in which R represents the cumulative rainfall.The second common method is Pruned Sequential Time Series (see for example [24,25]), which retains only the most relevant period of data sequence as the model input to prevent excessive variables from influencing the model simulations.For example, R(t − a), . . ., R(t − b), where t − a through t − b represent the time period when R is most correlated to the water level.
In addition to the above two methods, a third method namely Non-Sequential Time Series has been proposed by Talei et al. [26].In their study, two antecedent rainfall data R(t − T 1 ) and R(t − T 2 ) are selected as the model input, in which T 1 and T 2 are two non-sequential time points determined through a search test.Their results show that the models acquired in this method perform better than those acquired in the two aforementioned methods.In the study of Talei et al. [26], T 1 and T 2 were determined by a search process using the method of exhaustion.However, this approach of exhaustion search requires tremendous computing resources and CPU time, which makes it only suitable for the combination of a few regressors.The present study improves this aspect by adopting genetic algorithm (GA) for the optimal selection of model regressors, such that the number of model regressors is no longer limited and a vast amount of models with various combinations of regressors can be explored.
In the present study, the input variables of the models are selected non-sequentially from the combination of R(t), R(t − 1), . . ., R(t − m), where R represents the 18-h cumulative rainfall which, according to the previous data analysis, has been shown most correlated to the water-level.m defines the range of regressors to be selected.The greater the value of m, the more possible combinations of regressors and thus candidate models can be explored.Considering the computing time required for the search process, m is set to be 10 in the present study.The water level forecast is also related to water level data H(t), . . ., H(t − n).For simplicity, n is also set to be 10 herein.The relation between the regressors and output of each model is summarized below: , and H(t), H(t − 1), • • • , H(t − 10)] (7) in which ψ represents the function of ARX.

Assessing Indexes
In order to search for the optimal models, the following indexes are employed to evaluate the water level forecast capacity of each candidate model.

Coefficient of Efficiency (CE)
CE is an index designed to evaluate the performance of a hydrological model [27].CE is defined as follows where H obs and H est are the observed and estimated water levels, respectively; H obs is the average of observed water level, and n t is the number of data points.A CE value closer to 1 indicates that the predicted water levels fit more for the observations.

Relative Time Shift Error (RTS)
It has been shown by Talei and Chua [28] that a certain time shift error might emerge when using past data to forecast into the future.The forecasted water-level hydrograph displays a shifted delay in time from the observations.In order to assess the time shift error of the models, the process of Talei and Chua [28] is adopted in this study.The water-level hydrograph predicted by the model is first shifted back in time from 0 to k data points, and the CE of each shifted hydrograph is then calculated by comparing to the observations.The shift step δ associated to the maximum CE is considered as the time shift error of the model.The definition of Relative Time Shift is thus as follows where δ denotes the average δ of the typhoon events, and k is the forecast horizon (i.e., prediction time steps; each time step is 10 min in the present study).RTS ranges between 0 and 1.A smaller RTS represents less time shift in the predictions.
Threshold Statistic for a Level of x% (TS x ) TS x ( [29,30]) was employed to assess the error distribution in the forecasted water-level hydrographs.TS x is defined as follows where n is the total amount of data points; y x is the amount of forecasted data points with absolute relative error |RE t | less than a specified criteria of x%.RE t is defined as where H t o and H t c are the observed and predicted water level at time t, respectively; and d t o is the observed water depth.Using water depth instead of water level as the denominator in Equation ( 11) has several benefits.First, the value of |RE t | thus defined ranges between 0 and 1, where |RE t | = 0 indicates perfect prediction and 1 denotes a bad prediction deviated from the observation by the entire water depth.Second, a better accuracy in error measurement is achieved since the scale of water level might be in tens or hundreds of meters (depends on the location where the data were recorded) whereas in water depth it is often in meters only.A greater TS x means that the forecasted water-level hydrograph contains more predicted data points with the absolute relative error less than x%, and thus the model performance is better.The present study adopts 15% as the threshold value (i.e., TS 15 ).

Data Processing
The water level and cumulative rainfall data are standardized using the following equation [31] y n = 0.1 + 0.8 y i − y min y max − y min (12) where y n is the data after standardization; y i is the raw data; y max and y min are the maximum and minimum of raw data, respectively.The process of data standardization concentrates the dispersive data in a defined interval, which has shown great improvement on the predictions [31].

Model Optimization
The three indexes of CE, RTS, and TS 15 each evaluates the model's water level forecast accuracy, the time shift error, and the error distribution, respectively.All the three indexes provide valuable information for disaster prevention action during typhoon attack and shall be considered at the same time.However, it is difficult to weigh the importance of the three indexes.In order to find the models with good performance in all aspects, the multi-objective genetic algorithm (MOGA) was adopted in the present study as a tool for the search.

Multi-Objective Genetic Algorithm
The theoretical basis of GA was developed based on the Darwinian natural selection theory.After Holland [32] developed a firm mathematical foundation for the algorithm, GA has been widely applied to various fields to solve optimization problems that could not be tackled by traditional methods.In GA, each individual in a population is deemed as a possible solution to the optimization problem at hand.Based on the specified objective functions, the performance of each individual is evaluated and compared, and the individuals with better performance shall have a greater chance to pass its gene to the next generation.Through this procedure, the overall performance of the entire population gradually evolves and improves.After several generations of evolution, the individuals with optimal performance shall emerge, and these optimal individuals are considered as the optimal solutions of the problem.The algorithm has been shown capable of locating the global optima [33] and is particularly suitable for solving multi-objective optimization problems [34].

Objective Functions
Based on the definition of CE, RTS, and TS 15 , three design goals of the optimal model can be identified, which are the maximum CE, minimum RTS, and maximum TS 15 , respectively.Also, according to Equations ( 8)-(10), the upper limit of CE and TS 15 is 1, and RTS is always a positive number; therefore, the objective functions shall be defined as follows: Objective Objective 2 : minimize RTS ( Objective 3 : minimize 1 − TS 15 (15) in which CE, RTS, TS 15 represent the averaged value of the three indexes for the typhoon events, respectively.The design goal of Objective 1 and Objective 3 is to get CE and TS 15 closest to 1, while the design goal of Objective 2 is to get RTS closest to 0.

Codification of Regressor Combinations
The design variables of the optimization problem are the various possible combinations of model regressors consisting of R(t) through R(t − 10) and H(t) through H(t − 10), as shown in Equation (7).To utilize MOGA to search for the optimal models, the combination of the regressors has to be codified into a chromosome.In the present study, this procedure is accomplished by using binary bit string codifications, as illustrated in Figure 4.The combination of regressors is represented by a binary bit string chromosome consisting of 22 genes each associated to a specific regressor.If the value of a gene is 1, the regressor associated to this gene shall be selected as an input variable of the model.Contrarily, for a gene with a value of 0, the associated regressor shall not be selected as input.Each chromosome represents a specific combination of regressors and thus a candidate model.With 22 genes in one chromosome, the total candidate models thus formed in the search space are 2 22 ≈ 4.19 × 10 6 .
Water 2017, 9, 519 8 of 15 evaluated and compared, and the individuals with better performance shall have a greater chance to pass its gene to the next generation.Through this procedure, the overall performance of the entire population gradually evolves and improves.After several generations of evolution, the individuals with optimal performance shall emerge, and these optimal individuals are considered as the optimal solutions of the problem.The algorithm has been shown capable of locating the global optima [33] and is particularly suitable for solving multi-objective optimization problems [34].

Objective Functions
Based on the definition of CE, RTS, and TS 15 , three design goals of the optimal model can be identified, which are the maximum CE, minimum RTS, and maximum TS 15 , respectively.Also, according to Equations ( 8)-(10), the upper limit of CE and TS 15 is 1, and RTS is always a positive number; therefore, the objective functions shall be defined as follows: Objective in which CE , RTS , TS represent the averaged value of the three indexes for the typhoon events, respectively.The design goal of Objective 1 and Objective 3 is to get CE and TS closest to 1, while the design goal of Objective 2 is to get RTS closest to 0.

Codification of Regressor Combinations
The design variables of the optimization problem are the various possible combinations of model regressors consisting of ( ) through ( − 10) and ( ) through ( − 10), as shown in Equation (7).To utilize MOGA to search for the optimal models, the combination of the regressors has to be codified into a chromosome.In the present study, this procedure is accomplished by using binary bit string codifications, as illustrated in Figure 4.The combination of regressors is represented by a binary bit string chromosome consisting of 22 genes each associated to a specific regressor.If the value of a gene is 1, the regressor associated to this gene shall be selected as an input variable of the model.Contrarily, for a gene with a value of 0, the associated regressor shall not be selected as input.Each chromosome represents a specific combination of regressors and thus a candidate model.With 22 genes in one chromosome, the total candidate models thus formed in the search space are 2 4.19 × 10 .For model calibration, 6 out of the 10 events were selected, and the rest events were employed for validation.The selection of the calibration events is based on total cumulative rainfall to represent large, middle, and small amounts of rainfalls.The events thus selected are Songda, Nanmadol, Soulik, Fung-wong, Soudelor, and Dujuan.The three performance indexes of CE, RTS and TS 15 are For model calibration, 6 out of the 10 events were selected, and the rest events were employed for validation.The selection of the calibration events is based on total cumulative rainfall to represent large, middle, and small amounts of rainfalls.The events thus selected are Songda, Nanmadol, Soulik, Fung-wong, Soudelor, and Dujuan.The three performance indexes of CE, RTS and TS 15 are calculated and averaged over the validation events.The optimal models for the three types of LARX, NLARX-W, and NLARX-S are searched for using MOGA.With regard to MOGA setting, the population size is set to 200 with Pareto fraction of 0.35.The maximum generation of evolution is 500, and the stopping criterion of the evolution is set to 50 generations of stalls.

Results and Discussion
The result acquired by MOGA is the Pareto optimal model set, in which every model is un-dominated, that is, at least one of the three indexes of the model is not exceeded by another model.For each model type of LARX, NLARX-W, and NLARX-S, the three models with the best performance in the three indexes, respectively, are selected among the Pareto optimal set.The results are shown in Table 2.The selected models are named according to their model type and the featuring index, L1 represents the model with the maximum CE in LARX, L2 represents the one with the minimum RTS in LARX, and L3 represents the maximum TS 15 model in LARX.Similarly, W1, W2 and W3, as well as S1, S2, and S3, respectively represent the models with the maximum CE, minimum RTS, and maximum TS 15 in NLARX-W and NLARX-S.The scores of the nine optimal models on the three indexes and the corresponding chromosomes (i.e., selected regressors) are listed in Table 2.As shown, the selected regressors associated to each of the nine optimal models appear to be in a non-sequential pattern.This supports the result of Talei et al. [26] that a model with non-sequential inputs has better performance than the ones with sequential or pruned sequential inputs.The prediction lead time is crucial for disaster relief action during typhoon attack.An appropriate choice of the prediction lead time has to be linked to the properties of the hydrological behavior of the watershed.It has been observed in the study area that the variation of water level often lags behind the rainfall (for example, as observed in Figure 2).This time lag behavior between the rainfall and the water level presents a characteristic of the watershed, which can be used as an index for the selection of the appropriate prediction lead time.For that, the corrections between the cumulative rainfall data and the water level data shifted back in time by various time lags are analyzed.The results are as shown in Figure 5.The circle dots represent the average CC of the typhoon events and the error bars denote the maximum and the minimum CCs of the events.As seen, the average CC reaches the peak at the time lag of t − 3.This indicates that the water level is most correlated to the cumulative rainfall with 3 h of lead in the study area.The prediction lead time in the present study is thus selected to be 3 h in accordance with this hydrological characteristic of the watershed.In the comparison of CE, as shown in Figure 6a, the performance of Model Type 1 (L1, W1, S1) and Model Type 3 (L3, W3, S3) is quite good, in which CE all reaches above 0.8.The CEs of Type 1 Models are seen a little higher than Type 3 Models, but the differences are not significant.The CE performance of Model Type 2 is relatively poorer, especially for L2 and W2.This is because Model Type 2 features on reducing the time shift error.Nevertheless, it is noted that the CE of Model S2 still reaches above 0.7.As shown in Figure 6a, the CEs of the nonlinear models (W-and S-series of models) are somewhat higher than the linear L-series of models, thus indicating that the relationship between rainfall and water level at the site of WG2 is nonlinear.In the comparison of nonlinear models, the CEs of NLARX-S models (S1, S2, S3) appear to be higher than NLARX-W models (W1, W2, W3), and that difference is particularly evident with S2 and W2.In the comparison of time shift errors, as shown in Figure 6b, the RTS of Model Type 2 (L2, W2, S2) is seen lower than Type 1 (L1, W1, S1) and Type 3 (L3, W3, S3).The RTS of S2 in all models is the lowest, showing the minimum time shift error in the forecast, which is followed by W2 with a little increase in RTS; both of which are nonlinear models.The RTS of L2 is the worst in Type 2 models, and as shown in Figure 6b, even W1 and S1 which feature on CE have shown somewhat better performance on RTS than L2.Also, the RTS of the linear L1 and L3 models are much higher than the nonlinear W-and S-series of models.This yet again implies the nonlinearity between the rainfall and the water-level at the study area.Figure 6c is the comparison of the nine models in TS 15 performance.As shown in the figure, both Type 1 models ((L1, W1, S1) and Type 3 models (L3, W3, S3) display high performance on TS 15 .Among the nine models, the TS 15 of S3 is the highest, which is closely followed by W3.The linear L3 has the worst TS 15 in Type 3 models, which as shown in Figure 6c, is even lower than the nonlinear W1 and S1 that features on CE.In summary, the comparisons in Figure 6 show that the overall performance of nonlinear models is better than linear models on every aspect.Comparisons in the nonlinear models show that the NLARX-S models perform a little better than the NLARX-W models, but the differences are not very significant.
The results of the model predictions might be related to the hydrological behavior of the watershed.As has been shown in Figure 5, the correction between the water level and the rainfall in the study area is most prominent with a time lag of 3 h.It is noted in Table 2 that the model regressors optimized by MOGA all include R(t − 3) in the inputs.This result might reflect the hydrological behavior of the study area.It is also noted that Type 1 and Type 3 models which exhibit higher CE and TS15 values include not only R(t − 3) but also R(t − 4) or R(t − 2) as their regressors.As seen in  ) is quite good, in which CE all reaches above 0.8.The CEs of Type 1 Models are seen a little higher than Type 3 Models, but the differences are not significant.The CE performance of Model Type 2 is relatively poorer, especially for L2 and W2.This is because Model Type 2 features on reducing the time shift error.Nevertheless, it is noted that the CE of Model S2 still reaches above 0.7.As shown in Figure 6a, the CEs of the nonlinear models (W-and S-series of models) are somewhat higher than the linear L-series of models, thus indicating that the relationship between rainfall and water level at the site of WG2 is nonlinear.In the comparison of nonlinear models, the CEs of NLARX-S models (S1, S2, S3) appear to be higher than NLARX-W models (W1, W2, W3), and that difference is particularly evident with S2 and W2.In the comparison of time shift errors, as shown in Figure 6b, the RTS of Model Type 2 (L2, W2, S2) is seen lower than Type 1 (L1, W1, S1) and Type 3 (L3, W3, S3).The RTS of S2 in all models is the lowest, showing the minimum time shift error in the forecast, which is followed by W2 with a little increase in RTS; both of which are nonlinear models.The RTS of L2 is the worst in Type 2 models, and as shown in Figure 6b, even W1 and S1 which feature on CE have shown somewhat better performance on RTS than L2.Also, the RTS of the linear L1 and L3 models are much higher than the nonlinear W-and S-series of models.This yet again implies the nonlinearity between the rainfall and the water-level at the study area.Figure 6c is the comparison of the nine models in TS 15 performance.As shown in the figure, both Type 1 models ((L1, W1, S1) and Type 3 models (L3, W3, S3) display high performance on TS 15 .Among the nine models, the TS 15 of S3 is the highest, which is closely followed by W3.The linear L3 has the worst TS 15 in Type 3 models, which as shown in Figure 6c, is even lower than the nonlinear W1 and S1 that features on CE.In summary, the comparisons in Figure 6 show that the overall performance of nonlinear models is better than linear models on every aspect.Comparisons in the nonlinear models show that the NLARX-S models perform a little better than the NLARX-W models, but the differences are not very significant.The results of the model predictions might be related to the hydrological behavior of the watershed.As has been shown in Figure 5, the correction between the water level and the rainfall in the study area is most prominent with a time lag of 3 h.It is noted in Table 2 that the model regressors optimized by MOGA all include R(t − 3) in the inputs.This result might reflect the hydrological behavior of the study area.It is also noted that Type 1 and Type 3 models which exhibit higher CE and TS 15 values include not only R(t − 3) but also R(t − 4) or R(t − 2) as their regressors.As seen in Figure 2, the correction between the water level and the rainfall in the study area is also high with time lags of t − 2 and t − 4.This might explain the rather good results of CE and TS 15 scores achieved by Type 1 and Type 3 models.Figure 7 is the comparison of validated water level hydrograph and measured data of Type 1 models, with a prediction lead time of 3 h.Figure 7a-d show the validation results of Typhoon Saola, Matmo, Trami and Usagi, respectively.As seen in the figures, the models generally give reasonable forecasts comparing to the data.It is noted that all the three models exhibit certain degrees of time shift errors on the rising limb of the hydrographs.Since the rising limb of the flood wave is the most important phase for flood forecasting activities, the delayed forecasts pose a limitation of the models.Figure 7 is the comparison of validated water level hydrograph and measured data of Type 1 models, with a prediction lead time of 3 h.Figure 7a-d show the validation results of Typhoon Saola, Matmo, Trami and Usagi, respectively.As seen in the figures, the models generally give reasonable forecasts comparing to the data.It is noted that all the three models exhibit certain degrees of time shift errors on the rising limb of the hydrographs.Since the rising limb of the flood wave is the most important phase for flood forecasting activities, the delayed forecasts pose a limitation of the models.
Figure 7 is the comparison of validated water level hydrograph and measured data of Type 1 models, with a prediction lead time of 3 h.Figure 7a-d show the validation results of Typhoon Saola, Matmo, Trami and Usagi, respectively.As seen in the figures, the models generally give reasonable forecasts comparing to the data.It is noted that all the three models exhibit certain degrees of time shift errors on the rising limb of the hydrographs.Since the rising limb of the flood wave is the most important phase for flood forecasting activities, the delayed forecasts pose a limitation of the models.The above model performance comparison is mainly based on 3 h of prediction lead time.In order to investigate the performance under other lead times, the models are applied to the forecasts with prediction lead times varying from 0.5 h to 3 h, and the three indexes of CE, RTS, and TS 15 associated to each prediction lead time are calculated.The results are shown in Figure 8a-c.Figure 8a shows the CE variations of Type 1 models (L1, W1, and S1) along with the prediction lead.As shown in the figure, it appears that the CEs of all the three models increase with smaller prediction leads.For 0.5 h of prediction lead, the CEs of all models reach above 0.95, among which S1 is the highest.As the prediction lead time increases, the CEs of the three models gradually decrease.CE of L1 drops below 0.9 after the lead time passes 2 h, while W1 and S1 maintain above 0.9 up to 3 h of prediction lead.As the figure shows, for prediction leads from 0.5 to 3 h, CE performance of the three models shows S1 as optimal, closely followed by W1, and L1 is the worst.Figure 8b shows the RTS of Type 2 models (L2, W2, S2) varying along with prediction lead time from 0.5 to 3 h.As shown in the figure, the RTS of the three models gradually increases as the prediction lead time increases.All three models appear to have higher RTS with longer prediction leads, indicating greater relative time shift errors there.Comparing the three models, the RTS performance of S2 appears to be the best, and L2 is the worst.Figure 8c shows the TS 15 of Type 3 models (L3, W3, S3) varying along with prediction lead time.As shown in the figure, the TS 15 of all three models increase with smaller prediction leads.With 0.5 h of lead time, TS 15 of the three models all reach above 0.9, showing very good forecasts.As the prediction lead time increases, the TS 15 of the three models gradually decreases.The TS 15 performance of the three models shows S3 as the best, closely followed by W3, and L3 as the worst.The above model performance comparison is mainly based on 3 h of prediction lead time.In order to investigate the performance under other lead times, the models are applied to the forecasts with prediction lead times varying from 0.5 h to 3 h, and the three indexes of CE, RTS, and TS 15 associated to each prediction lead time are calculated.The results are shown in Figure 8a-c.Figure 8a shows the CE variations of Type 1 models (L1, W1, and S1) along with the prediction lead.As shown in the figure, it appears that the CEs of all the three models increase with smaller prediction leads.For 0.5 h of prediction lead, the CEs of all models reach above 0.95, among which S1 is the highest.As the prediction lead time increases, the CEs of the three models gradually decrease.CE of L1 drops below 0.9 after the lead time passes 2 h, while W1 and S1 maintain above 0.9 up to 3 h of prediction lead.As the figure shows, for prediction leads from 0.5 to 3 h, CE performance of the three models shows S1 as optimal, closely followed by W1, and L1 is the worst.Figure 8b shows the RTS of Type 2 models (L2, W2, S2) varying along with prediction lead time from 0.5 to 3 h.As shown in the figure, the RTS of the three models gradually increases as the prediction lead time increases.All three models appear to have higher RTS with longer prediction leads, indicating greater relative time shift errors there.Comparing the three models, the RTS performance of S2 appears to be the best, and L2 is the worst.Figure 8c shows the TS 15 of Type 3 models (L3, W3, S3) varying along with prediction lead time.As shown in the figure, the TS 15 of all three models increase with smaller prediction leads.With 0.5 h of lead time, TS 15 of the three models all reach above 0.9, showing very good forecasts.As the prediction lead time increases, the TS 15 of the three models gradually decreases.The TS 15 performance of the three models shows S3 as the best, closely followed by W3, and L3 as the worst.

Date
shift errors there.Comparing the three models, the RTS performance of S2 appears to be the best, and L2 is the worst.Figure 8c shows the TS 15 of Type 3 models (L3, W3, S3) varying along with prediction lead time.As shown in the figure, the TS 15 of all three models increase with smaller prediction leads.With 0.5 h of lead time, TS 15 of the three models all reach above 0.9, showing very good forecasts.As the prediction lead time increases, the TS 15 of the three models gradually decreases.The TS 15 performance of the three models shows S3 as the best, closely followed by W3, and L3 as the worst.

Conclusions
A methodology for the determination of the optimal combination of non-sequential regressors for typhoon inundation forecasting models has been developed.By integrating a MOGA with ARXbased models, the proposed methodology is capable of locating the optimal models that conform to multiple objectives in terms of high prediction accuracy, low time shift error, and low threshold statistics.Testing results show that the nine optimal models obtained by the MOGA all display nonsequential patterns in the resultant combinations of regressors, which signifies the superiority of this type of model as well as the capacity of the proposed methodology in locating the optimal nonsequential regressors.In comparing the resultant models obtained by the proposed methodology, the results show that the overall performance of nonlinear models (NLARX-W and NLARX-S) is significantly better than linear models (LARX), revealing a nonlinear relationship between the rainfall and water level at the study area.On all the three assessing indexes of CE, RTS, and TS 15 , the evaluation results of the three types of models present, in general, NLARX-S as the best, NLARX-W falling slightly behind, and LARX as the worst.These results provide a new approach to constructing better models for typhoon inundation level forecasts.

Conclusions
A methodology for the determination of the optimal combination of non-sequential regressors for typhoon inundation forecasting models has been developed.By integrating a MOGA with ARX-based models, the proposed methodology is capable of locating the optimal models that conform to multiple objectives in terms of high prediction accuracy, low time shift error, and low threshold statistics.Testing results show that the nine optimal models obtained by the MOGA all display non-sequential patterns in the resultant combinations of regressors, which signifies the superiority of this type of model as well as the capacity of the proposed methodology in locating the optimal non-sequential regressors.In comparing the resultant models obtained by the proposed methodology, the results show that the overall performance of nonlinear models (NLARX-W and NLARX-S) is significantly better than linear models (LARX), revealing a nonlinear relationship between the rainfall and water level at the study area.On all the three assessing indexes of CE, RTS, and TS 15 , the evaluation results of the three types of models present, in general, NLARX-S as the best, NLARX-W falling slightly behind, and LARX as the worst.These results provide a new approach to constructing better models for typhoon inundation level forecasts.

Figure 2 .
Figure 2. Water level and rainfall records during Typhoon Trami in Donsan area, Taiwan.

Figure 2 .
Figure 2. Water level and rainfall records during Typhoon Trami in Donsan area, Taiwan.
Water 2017, 9, 519 5 of 15Two different types of nonlinear ARX models are employed in the present study, viz. the nonlinear ARX with Wavelet function (NLARX-W) and the nonlinear ARX with Sigmoid function (NLARX-S).In NLARX-W, the nonlinear unit adopts the Wavelet function[17], which has the following formula:

Figure 3 .
Figure 3. Correlations between water level and cumulative rainfall with various accumulation durations.

Figure 3 .
Figure 3. Correlations between water level and cumulative rainfall with various accumulation durations.

Figure 4 .
Figure 4. Binary bit string codification of a combination of regressors.

Figure 4 .
Figure 4. Binary bit string codification of a combination of regressors.

Water 2017, 9 , 519 10 of 15 lag of t − 3 .
This indicates that the water level is most correlated to the cumulative rainfall with 3 h of lead in the study area.The prediction lead time in the present study is thus selected to be 3 h in accordance with this hydrological characteristic of the watershed.

Figure 5 .
Figure 5. Correlations between cumulative rainfall and water level with various time lags.

Figure
Figure6a-c compare the performance of the nine models in CE, RTS, and TS 15 .In the comparison of CE, as shown in Figure6a, the performance of Model Type 1 (L1, W1, S1) and Model Type 3 (L3, W3, S3) is quite good, in which CE all reaches above 0.8.The CEs of Type 1 Models are seen a little higher than Type 3 Models, but the differences are not significant.The CE performance of Model Type 2 is relatively poorer, especially for L2 and W2.This is because Model Type 2 features on reducing the time shift error.Nevertheless, it is noted that the CE of Model S2 still reaches above 0.7.As shown in Figure6a, the CEs of the nonlinear models (W-and S-series of models) are somewhat higher than the linear L-series of models, thus indicating that the relationship between rainfall and water level at the site of WG2 is nonlinear.In the comparison of nonlinear models, the CEs of NLARX-S models (S1, S2, S3) appear to be higher than NLARX-W models (W1, W2, W3), and that difference is particularly evident with S2 and W2.In the comparison of time shift errors, as shown in Figure6b, the RTS of Model Type 2 (L2, W2, S2) is seen lower than Type 1 (L1, W1, S1) and Type 3 (L3, W3, S3).The RTS of S2 in all models is the lowest, showing the minimum time shift error in the forecast, which is followed by W2 with a little increase in RTS; both of which are nonlinear models.The RTS of L2 is the worst in Type 2 models, and as shown in Figure6b, even W1 and S1 which feature on CE have shown somewhat better performance on RTS than L2.Also, the RTS of the linear L1 and L3 models are much higher than the nonlinear W-and S-series of models.This yet again implies the nonlinearity between the rainfall and the water-level at the study area.Figure6cis the comparison of the nine models in TS 15 performance.As shown in the figure, both Type 1 models ((L1, W1, S1) and Type 3 models (L3, W3, S3) display high performance on TS 15 .Among the nine models, the TS 15 of S3 is the highest, which is closely followed by W3.The linear L3 has the worst TS 15 in Type 3 models, which as shown in Figure6c, is even lower than the nonlinear W1 and S1 that features on CE.In summary, the comparisons in Figure6show that the overall performance of nonlinear models is better than linear models on every aspect.Comparisons in the nonlinear models show that the NLARX-S models perform a little better than the NLARX-W models, but the differences are not very significant.The results of the model predictions might be related to the hydrological behavior of the watershed.As has been shown in Figure5, the correction between the water level and the rainfall in the study area is most prominent with a time lag of 3 h.It is noted in Table2that the model regressors optimized by MOGA all include R(t − 3) in the inputs.This result might reflect the hydrological behavior of the study area.It is also noted that Type 1 and Type 3 models which exhibit higher CE and TS15 values include not only R(t − 3) but also R(t − 4) or R(t − 2) as their regressors.As seen in

Figure 5 .
Figure 5. Correlations between cumulative rainfall and water level with various time lags.

Figure
Figure6a-c compare the performance of the nine models in CE, RTS, and TS 15 .In the comparison of CE, as shown in Figure6a, the performance of Model Type 1 (L1, W1, S1) and Model Type 3 (L3, W3, S3) is quite good, in which CE all reaches above 0.8.The CEs of Type 1 Models are seen a little higher than Type 3 Models, but the differences are not significant.The CE performance of Model Type 2 is relatively poorer, especially for L2 and W2.This is because Model Type 2 features on reducing the time shift error.Nevertheless, it is noted that the CE of Model S2 still reaches above 0.7.As shown in Figure6a, the CEs of the nonlinear models (W-and S-series of models) are somewhat higher than the linear L-series of models, thus indicating that the relationship between rainfall and water level at the site of WG2 is nonlinear.In the comparison of nonlinear models, the CEs of NLARX-S models (S1, S2, S3) appear to be higher than NLARX-W models (W1, W2, W3), and that difference is particularly evident with S2 and W2.In the comparison of time shift errors, as shown in Figure6b, the RTS of Model Type 2 (L2, W2, S2) is seen lower than Type 1 (L1, W1, S1) and Type 3 (L3, W3, S3).The RTS of S2 in all models is the lowest, showing the minimum time shift error in the forecast, which is followed by W2 with a little increase in RTS; both of which are nonlinear models.The RTS of L2 is the worst in Type 2 models, and as shown in Figure6b, even W1 and S1 which feature on CE have shown somewhat better performance on RTS than L2.Also, the RTS of the linear L1 and L3 models are much higher than the nonlinear W-and S-series of models.This yet again implies the nonlinearity between the rainfall and the water-level at the study area.Figure6cis the comparison of the nine models in TS 15 performance.As shown in the figure, both Type 1 models ((L1, W1, S1) and Type 3 models (L3, W3, S3) display high performance on TS 15 .Among the nine models, the TS 15 of S3 is the highest, which is closely followed by W3.The linear L3 has the worst TS 15 in Type 3 models, which as shown in Figure6c, is even lower than the nonlinear W1 and S1 that features on CE.In summary, the comparisons in Figure6show that the overall performance of nonlinear models is better than linear models on every aspect.Comparisons in the nonlinear models show that the NLARX-S models perform a little better than the NLARX-W models, but the differences are not very significant.

Water 2017, 9 , 519 11 of 15 Figure 6 .
Figure 6.Performance comparison of the optimal models (a) Coefficient of efficiency (CE), (b) Relative time shift error (RTS) and (c) Threshold statistic for a level of % TS 15 .

Figure 6 .
Figure 6.Performance comparison of the optimal models (a) Coefficient of efficiency (CE), (b) Relative time shift error (RTS) and (c) Threshold statistic for a level of x% TS 15 .

Figure 8 .
Figure 8. Variation of assessing indexes of the optimal models with respect to different prediction lead times (a) CE, (b) RTS and (c) TS15.

Figure 8 .
Figure 8. Variation of assessing indexes of the optimal models with respect to different prediction lead times (a) CE, (b) RTS and (c) TS 15 .

Table 1 .
Historical typhoon records in Donsan area.