JP7749341B2 - Modular time series data prediction device, modular time series data prediction method, and program - Google Patents

Description

Embodiments of the present invention relate to a modular time series data prediction device, a modular time series data prediction method, and a program.

Infrastructure systems such as water supply and sewerage systems, stormwater drainage systems, power systems, and transportation systems, as well as industrial process plants such as steel processes, petrochemical processes, and semiconductor manufacturing processes, typically have multiple online sensors installed to measure multiple process conditions. A system known as a process monitoring and control system (SCADA: Supervisory Control And Data Acquisition) acquires process data (flow rate, temperature, water quality, manipulated variables, etc.) obtained through measurements by sensors installed in the above infrastructure systems and industrial plants, and stores it on a server as time-series data. Process monitoring and control systems typically provide this time-series data to supervisors as trend graphs.

In addition to directly monitoring trends in time-series data as described above, monitors may also monitor plants using predictive information. Plant monitoring systems that use predictive information can predict future values of the monitored object, in addition to real-time actual measurements (measured values) during monitoring, and provide the predicted results to operation managers, thereby helping them to implement appropriate operation plans. Plant monitoring systems that use predictive information are a typical means of advanced monitoring, along with abnormal sign detection and anomaly diagnosis in plants.

Examples of infrastructure systems that use plant monitoring that utilizes predictive information include predicting water demand in water purification systems, predicting wastewater inflow into sewage treatment plants, predicting rainwater inflow into sewage treatment plants and stormwater drainage pumping stations during rainy weather, predicting river and dam water levels during rainy weather, short-term rainfall forecasts (nowcasting) using weather radar, traffic volume forecasts and congestion forecasts in transportation systems, power demand forecasts in power systems, and customer numbers for cafeterias and shops within facilities such as buildings. These systems are widely used in a variety of infrastructure fields.

Plant monitoring systems that use predictive information not only provide the predictive information itself as added value, but can also be used to achieve more efficient and safer operations by using the predictive information in plant operation plans and plant control. For example, demand forecasts for water demand, electricity demand, etc. can be used to optimize water operation plans and power generation and storage plans (EMS: Energy Management Systems) to achieve more efficient operations. Furthermore, for example, rainfall forecasts and inflow forecasts can be used as input information for controlling the operation of stormwater drainage pumps, contributing to risk reduction by preventing or suppressing flooding.

Time series data analysis is often used not only in infrastructure systems but also in process plants such as steel and petrochemical processes to predict product quality and yield, and this information is often used to improve production efficiency (yield).

Time-series data forecasting plays an extremely important role in infrastructure systems and process industrial plants, and various data-driven (black-box approach) forecasting methods have traditionally been applied using plant monitoring data.

Prediction techniques using time-series data such as plant monitoring data have been widely used for some time. However, with the rapid advances in AI-related theories, technologies, and methods in recent years, particularly in the field of statistical machine learning, AI methods for time-series data analysis have also been advancing in parallel, accelerating the trend toward applying various advanced AI methods to time-series data prediction, regardless of field or area. In particular, deep learning networks (DNNs), an advanced form of neural networks, are a representative method that has been attracting attention in recent years, and DNNs are sometimes used as a synonym for AI.

In addition to DNN, a representative AI method, there are various methods based on regularization that were introduced (from a different perspective than the previously mentioned PLS) as a way to avoid instability (ill-conditioned problem) caused by multicollinearity, which is a problem in multiple regression analysis. In particular, in addition to the long-standing L2 regularization (ridge regression), there are methods such as Lasso, an L1 regularization method that can automatically select the necessary explanatory variables from a large number of candidate explanatory variables, and RVR (relevant vector regression), a regression method (prediction method) based on RVM (relevant vector machine) that automatically selects explanatory variables using a Bayesian inference method from a different perspective than Lasso, and methods that use optimization to solve classification problems ( A variety of advanced time series data forecasting methods are also coming into use, such as SVR (support vector regression), which adapts SVMs (support vector machines), which have been made robust for classification problems, to regression problems (prediction problems); Gaussian process regression, which nonlinearly extends linear regression such as conventional multiple regression and further considers probability distributions to reduce it to a nonparametric regression problem, outputting a predictive distribution rather than a point prediction; and various methods known as bagging and boosting (random forest, adaboost, XG boost, etc.), which combine various prediction subsystems to improve predictive performance.

In this way, advanced AI-based prediction methods are being developed one after another, and there have been many reports of significantly improved prediction accuracy compared to traditional methods, raising expectations for AI-related technologies.

On the other hand, one of the reasons why AI-related technologies are not necessarily easily deployed in actual, realistic systems (infrastructure systems or process plants) is the problem of "lack of explainability of results," which has been pointed out as a drawback of the black-box approach that is an essential characteristic of statistical and AI analysis methods. Awareness of this problem of "lack of explainability of results" is rapidly increasing, and in an effort to solve this problem, the concept of "explainable AI (XAI: eXplainable AI)" has recently been proposed by the US DARPA, which has attracted attention.

Patent No. 6261960

Infrastructure systems often use white-box models based on physical and chemical knowledge, which is thought to be related to the need for "accountability." In fact, the predictive accuracy of white-box models is not necessarily higher than that of black-box models, and white-box models do not necessarily accurately and faithfully reproduce actual physical phenomena. Despite this, white-box models are sometimes preferred and used because the "parameters" contained in white-box models have physical and chemical meanings, making them "easy to explain rationally and easy to convince others." Furthermore, because these parameters have physical and chemical meanings, it is easy to establish policies and guidelines for how to adjust them if accuracy is insufficient. In fact, parameters with physical and chemical meanings are easy to understand because the parameter values themselves have physical and chemical meanings. It is possible to imagine the range of possible values for the parameter and consider the impact of increasing or decreasing the value on the results based on physical and chemical knowledge.

As mentioned above, "ease of rational explanation and ease of adjustment" may be a major motivation for using white box models. If the elements of "ease of rational explanation and ease of adjustment" could be added to black box models using AI-like techniques, it is believed that their range of applications would be greatly expanded.

Embodiments of the present invention have been made in consideration of the above circumstances, and aim to provide a modular time series data prediction device, a modular time series data prediction method, and a program that rationally explain prediction results and allow for easy adjustment of prediction results.

A modular time-series data prediction device according to an embodiment is a device applied to a system or process having a plurality of process sensors that measure a plurality of process variables at a predetermined cycle, and includes a data collection and storage unit that collects and stores time-series data of the plurality of process variables at a predetermined cycle, and that includes an output variable data selection unit that selects at least one output variable to be predicted from the plurality of process variables, and an input variable data selection unit that selects a plurality of input variables from the plurality of process variables; and a data collection and storage unit that performs pairwise prediction for each one-input-one-output pair using identification data of the output variable and the plurality of input variables extracted from the time-series data. a pairwise prediction model identification unit that identifies parameters of a prediction model to define a plurality of pairwise prediction models; a prediction model synthesis method definition unit that defines a synthesis method for pairwise prediction values output from the plurality of pairwise prediction models; a pairwise output variable prediction unit that inputs prediction data for the plurality of input variables extracted from the time-series data in a predetermined cycle in the time progression direction or in real time to the plurality of pairwise prediction models to calculate the pairwise prediction values corresponding to each of the plurality of input variables; and a synthesis output variable prediction unit that synthesizes the plurality of pairwise prediction values by the synthesis method to calculate a prediction value of the output variable.

FIG. 1 is a diagram illustrating a rainfall inflow prediction system to which a modular time-series data prediction device according to an embodiment of the present invention is applied. FIG. 2 is a diagram illustrating an example of the configuration of a modular time-series data prediction device according to an embodiment. FIG. 3 is a scatter diagram that schematically illustrates an example of the relationship between main water level data and inflow data. FIG. 4 is a scatter diagram that schematically illustrates an example of the relationship between main water level data and inflow data. FIG. 5 is a diagram illustrating a first example of a pairwise prediction model in a modular time series prediction device according to an embodiment. FIG. 6 is a diagram illustrating a second example of a pairwise prediction model in a modular time series prediction device according to an embodiment. FIG. 7 is a diagram for explaining the effect of the pairwise prediction model of the second embodiment. FIG. 8 is a diagram for explaining the effect of the pairwise prediction model of the second embodiment. FIG. 9 is a diagram illustrating a third example of a pairwise prediction model in a modular time series prediction device according to an embodiment. FIG. 10 is a diagram for explaining an example of the effect of the pairwise prediction model of the third embodiment. FIG. 11 is a diagram for explaining an example of the effect of the pairwise prediction model of the third embodiment.

In this embodiment, a modular time series data prediction device is described as a basic means for performing adjustable time series prediction. The modular time series prediction device of this embodiment allows a prediction model (system) to be decomposed into sub-models (subsystems), making each sub-model easily attachable and detachable.

A modular time series data prediction device can prevent the effects of, for example, when data for one of multiple monitored variables is missing, when multiple monitored variables contain outliers (abnormal data), or when monitored variables are unstable measurement data containing a large amount of noise, from spreading to the entire system. Furthermore, a modular time series data prediction device can easily determine whether explanatory variables are necessary (or not), and can easily detach partial models or partially adjust the prediction model.

These advantages are achieved by the modular structure, which allows the model to be "decomposable." The advantage of using a modular structure for predictive models is not simply that it is decomposable; in addition, this modular structure makes it possible to build a rational and convincing model that does not contradict the laws of physics, and this also makes it easier to adjust.

The advantages of a modular approach to predictive models are explained below.
It is generally believed that black-box approaches are not explainable because they only provide external information, such as inputs and outputs, and their internal structure is invisible. In other words, black-box approaches not only hide the internal structure, but also make it impossible to uniquely determine the internal structure from the external input-output relationships alone. There are potentially multiple (or even countless) internal structures that give the same input-output relationships. This makes it difficult to interpret the internal structure, and in some cases can lead to irrational models that contradict the laws of physics. This issue is known as the model uniqueness problem in the field of system identification, which is part of the field of control theory. The concept of model uniqueness is, for example, the question of whether, for a given pair of input and output data, it is possible to uniquely determine a model that represents the input-output relationship. This is known as the model identifiability problem.

The identifiability problem arises when the structure of a model (the form of the mathematical formula that represents the model) is limited in advance, for example, when the model is limited to a linear transfer function model or a time series model, and the question is whether the parameters included in that model (for example, the coefficients included in the transfer function) can be uniquely determined. In this case, the problem is specifically called a parameter identifiability problem. This parameter identifiability problem frequently arises in problems that deal with real-world data, and the following three cases are particularly commonly encountered.

The first is the problem of multicollinearity mentioned earlier, which occurs when there is a strong correlation between input variables (explanatory variables). This is a problem that typically arises when, for example, predicting the amount of rainwater inflow into a stormwater pumping station, using radar rainfall information as an input variable, and using the rainfall amounts of two different adjacent meshes as input variables; in real-world data, it is extremely common for there to be a strong correlation between explanatory variables.

The second problem occurs when the persistent exciting condition (PE condition), well known in the field of control theory, is not met. This occurs when input data that does not contain sufficient frequency components is applied to the process being modeled. An extreme, yet realistic, example is when the input data is constant, the output data is usually also constant. However, because there are countless transfer function models that represent this constant input-output relationship, it is not possible to uniquely determine the coefficients. For example, when considering a simple transfer function model (time series model) expressed as y(t) = b1 × u(t-1) + b2 × u(t-2), if u(t-1) = u(t-2) = k (a constant value), b1 and b2 cannot be uniquely determined. A generalization of this concept is the parameter identifiability condition known as the persistent exciting condition, but such conditions often do not hold.

The third case is when the amount of available data is small and the quality of the data is poor compared to the complexity of the model structure. A typical example is when it is known from a physical perspective that a certain object will exhibit nonlinear behavior, but the available data is insufficient, or only data is available around certain operating conditions, making it impossible to identify the nonlinear coefficients. For example, it is known that y(t) can be written as y(t)=c1×u(t)+c2×u(t)^2, but when the value of u(t) available as data is small, u(t)^2<<u(t), making it difficult to identify c2.

The above three examples are typical cases where parameter identifiability is lacking, but in real-world problems, such situations are extremely common. When modeling relatively small mechanical and electrical systems, it is common practice to conduct experiments to improve identifiability and build a model only after ensuring identifiability. However, when it comes to large-scale social infrastructure systems such as water supply and sewerage, stormwater drainage, and power systems, the data to be modeled is often data obtained from natural phenomena (rainfall, solar radiation, wind speed, etc.) or data from the target plant's operating state. This makes it impossible to conduct experiments to improve parameter identifiability, and the data available for modeling is often given, making it inherently difficult to avoid a lack of identifiability. As a result, when it comes to social infrastructure systems, models must be built without identifiability.

Hereinafter, the conventional technology and the modular time series data prediction device of this embodiment will be compared with regard to "multicollinearity," which is an example of a problem caused by a lack of identifiability.
"Multicollinearity" is a problem where the parameters of a regression model cannot be uniquely determined due to a correlation (linear dependence) between multiple explanatory variables (input variables). This will be explained using the simplest example.

Consider the following multiple regression model in which one output y is expressed by multiple regression of two inputs u1 and u2.
y = a1 × u1 + a2 × u2 (1)
Equation (1) is the simplest multiple regression model, with a1 and a2 being parameters called regression coefficients. To get a concrete physical image, let's consider the output y as the amount of rainwater inflow to a stormwater pumping station, and inputs u1 and u2 as the amounts of rainfall at different locations. In this case, it is self-evident from a physical perspective that an increase in rainfall will increase the amount of inflow, so it is easy to see that a1 and a2 must have positive values. Now, let's consider a case where there is a multicollinearity problem, where rainfall u1 and rainfall u2 are dependent on each other.

This multicollinearity problem is frequently observed in real-world phenomena. In this case, if we consider inputs u1 and u2 as rainfall amounts for adjacent meshes of radar rainfall data, it becomes clear that this problem is not unique, but rather a natural one. For example, in the weather radar called XRAIN provided by the Ministry of Land, Infrastructure, Transport and Tourism, the size of each mesh is 250m x 250m, so the rainfall amounts for adjacent meshes refer to rainfall amounts at locations several hundred meters apart, and it is quite natural from a physical perspective that these rainfall amounts are nearly identical. Even if the two rainfall data are not rainfall amounts for adjacent meshes of radar rainfall but ground rainfall data from two locations, it is easy to imagine that the values will be relatively close, given that the size of a stormwater drainage district is usually within a few square kilometers. Therefore, in such cases, it is natural that u1 ≒ u2, which is a true multicollinearity problem.

For simplicity's sake, let us assume that there is perfect multicollinearity, i.e., u1 = u2. For example, suppose that there is rainfall of 20 mm/h and rainwater inflow of 10 ^m3 /s. Now, since u1 = u2, if we define u: = u1 = u2, then equation (1) can be rewritten as equation (2) below.

y = (a1 + a2) × u (2)
The above equation (2) is in the form of a simple regression, and since u = u1 = u2, a1 + a2 = 0.5. However, as long as the constraint a1 + a2 = 0.5 is satisfied, equation (2) will produce the exact same output result regardless of the values of a1 and a2. Therefore, the values of a1 and a2 can have any value within the range that satisfies a1 + a2 = 0.5. Therefore, for example, it is no problem if a1 = 3.5 and a2 = -3 or if a1 = 100 and a2 = -99.5, and they give the exact same input-output relationship, so it becomes impossible to distinguish the differences in each parameter value from the input-output data.

On the other hand, since a1 and a2 are coefficients for each rainfall amount, they should essentially have at least positive values, and it is physically reasonable that if the rainfall amounts at two locations are the same, their values should be roughly equal. However, due to the correlation (linear dependence) between input variables, parameter values cannot be uniquely determined, and the true parameter values cannot be correctly identified; this is the essence of the problem of multicollinearity.

This lack of parameter identifiability, such as multicollinearity, makes it impossible to uniquely determine parameters from input-output relationships. In the fields of system identification and statistics (statistical machine learning = AI), this problem is treated as a mathematical "ill-posed problem," and various algorithms have been developed to avoid this ill-posedness. Methods such as the regularization and PLS mentioned above are typical technical means (techniques) for avoiding such ill-posedness, and many advanced AI methods employ techniques such as regularization and dimensionality reduction (used in PLS).

However, even if an advanced algorithm utilizing techniques such as regularization or PLS is applied to the problem in equation (1), it is not guaranteed that it will estimate values such that a1 > 0, a2 > 0, and a1 ≒ a2. In fact, it is highly likely that such values will not be obtained. This is because regularization and dimensionality reduction are techniques for avoiding the problem of multicollinearity (or its generalized equivalent, the lack of identifiability) as a mathematically ill-conditioned problem; they merely avoid ill-conditions and enable parameter estimation (identification). In other words, applying classical least squares to a multiple regression problem with multicollinearity results in the problem of "not being able to find an answer = not being able to find parameter values." Therefore, by introducing techniques for avoiding ill-conditions (regularization and dimensionality reduction) into the algorithm, the problem becomes "being able to find an answer = being able to identify parameter values," without directly considering whether the identified values are physically reasonable.

In fact, for example, in regularization techniques, in order to avoid bad conditions, when identifying parameters, not only is the prediction error evaluated, but an evaluation index is added that reduces the parameter norm (an index that measures magnitude) so that "the parameter value itself does not become too large." This is not directly related to assigning physical meaning or interpretation to the parameters.

On the other hand, a method for obtaining parameter values that are reasonably acceptable for such a problem will be described below. First, a simple regression model of equations (3) and (4) will be considered, instead of the model of equation (1).
y=a1´×u1 (3)
y=a2´×u2 (4)

In this case, when u1 = u2 and u1 = u2 = 20 mm/h, the model that outputs y = 10 ^m3 /s clearly has a1' = a2' = 0.5 if we consider the unit system as it is. Therefore, when such a relationship exists, if we apply simple regression to equations (3) and (4), a value close to a1' = a2' = 0.5 will be identified (learned) as the parameter. This value of 0.5 is a reasonable and convincing value (although the physical meaning is unclear because it is a regression model).

Next, consider constructing a multiple regression model of equation (1) using equations (3) and (4).
In this example, the input variables are the rainfall amounts at two locations, and since the sum of these is taken, the input rainfall amount doubles. Therefore, the coefficients a1 and a2 for each rainfall amount are half the coefficients of a simple regression model with one rainfall amount as input, and it is considered that a1 = a2 = 0.25 is a physically reasonable and plausible parameter value.

In this way, the relationship between u1 and y, and the relationship between u2 and y are determined using simple regression in equations (3) and (4), respectively, to construct two predictive models. Next, the average output y of the two predictive models is calculated. By following this constructive procedure, the values of a1 and a2 are each estimated to be 0.25, which are reasonable values that are considered physically reasonable, and as a result, the coefficients of the multiple regression equation in the form of equation (1) can be determined. This is the very essence of a modular approach, which combines partial models (subsystems).

The reason that the coefficients can be determined using the above multiple regression model is that, generally, when the same data is used to represent input/output relationships, models with simpler model structures (i.e., models with fewer parameters) have higher identifiability (i.e., parameter values can be uniquely identified). The modular time series data prediction device of this embodiment makes use of this characteristic. In other words, even when the same input/output data is used, identifiability depends on the model structure and number of parameters, and even if a model with a large number of parameters is not identifiable, it may be identifiable if the number of parameters is small.

In the above example, in the modular time series data forecasting device of this embodiment, a two-input multiple regression model becomes non-identifiable due to multicollinearity, but by considering this as a combination of one-input simple regressions, the problem of multicollinearity is essentially avoided, and rational model construction is made possible by stacking identifiable models.

Next, the relationship and differences between the modular time-series data prediction device of this embodiment and so-called "Explainable AI (XAI)" will be described.
So-called XAI is a top-down approach (or inductive approach) in which a black-box model is constructed by treating the whole as one large system, and then the model is broken down into an appropriately interpretable form and an explanation is added later.In contrast, the modular time-series data prediction device of this embodiment is a bottom-up approach (or deductive approach) in which an overall prediction model is constructed by integrating models (subsystems) that are individually relatively easy to understand and capable of rational explanation.

The bottom-up approach described above cannot take into account mutual interference between subsystems in advance, so it cannot be denied that it may be inferior to a top-down (XAI) approach when focusing on prediction accuracy. However, in terms of rational explainability, a bottom-up approach makes it possible to build a model that maintains parameter identifiability, as explained above, and is therefore likely to be superior to a top-down approach.

The bottom-up approach is also similar to the white-box modeling approach. White-box modeling is essentially a deductive, constructive modeling method. In white-box modeling, models are constructed constructively using parameters that have physical meaning, whereas in the modular approach of this embodiment, parameters that are easy to interpret are identified by directly relating each parameter to input and output, even though they cannot be given physical meaning, and models are constructed constructively using these parameters.

The following makes it clear that the modular approach of this embodiment is superior to the XAI approach in terms of rational explainability. For example, in the previous example, if a top-down (XAI) approach were to be used, the goal would be to establish a method for determining the coefficient parameters a1 and a2 using methods such as regularization or PLS, and then to add some kind of interpretation to the results to make them "explainable." However, as can be seen from the simple example above, if a1 and a2 are identified as values far removed from 0.25 (which is indeed possible when using regularization or PLS), it is easy to see that, in principle, a rational explanation is difficult to achieve, no matter what method is used to attempt an "explanation." For example, even if parameter values can be determined using PLS or regularization, if either a1 or a2 is a negative value (e.g., a1 = 3.5, a2 = -3), a rationally convincing explanation is impossible.

In other words, even if the value of 0.25, which is originally a reasonably acceptable value, was not obtained, in the absence of identifiability, as long as certain constraints are met, no difference will be seen as long as the focus is solely on prediction accuracy, regardless of the internal parameter values, and so the model will be considered to be of high accuracy. However, a high-accuracy model does not necessarily mean a rational and reasonable model that compensates for the validity of the model's parameter values, and so a top-down (XAI) approach makes it difficult in principle to provide a rational explanation.

In an XAI approach, it is difficult to assign meaning to a1 and a2 individually, but it is possible to assign meaning to a1 + a2, so it is not impossible to explain. However, in principle, it is highly likely that an XAI approach will be difficult to assign meaning to all of the parameters contained within, so some kind of "feature" will be extracted and a meaning will be sought for that.

Therefore, if the main goal is to build a "model with high rational explainability," it is highly likely that a bottom-up, constructive, deductive, modular approach that builds models to avoid states that lack identifiability and ensure the uniqueness of the model's internal structure (internal parameters), and integrates models while ensuring rational explainability, would be superior, rather than relying on algorithms that avoid mathematically ill-conditioning such as regularization or PLS.

In addition to being able to build a model that ensures rational explainability in this way, the separable structure mentioned at the beginning also makes it easier to adjust the model, and by combining these, it becomes possible to build a model that is rationally explainable and adjustable.
The modular time-series data prediction device of this embodiment is based on the above insight.

Hereinafter, a modular time-series data prediction device according to an embodiment will be described in detail with reference to the drawings.
FIG. 1 is a diagram illustrating a rainfall inflow prediction system to which a modular time-series data prediction device according to an embodiment of the present invention is applied.

The following describes a stormwater inflow prediction system to more clearly illustrate the implementation image and effects of a system equipped with the modular time series data prediction device of this embodiment. Note that the target process to which the modular time series data prediction device of this embodiment is applied is not limited to stormwater inflow prediction systems, but can be applied to any process in which time series data for two or more measurement items exists. The modular time series data prediction device can be applied to processes such as steel plants, petrochemical plants, food plants, pharmaceutical processes, semiconductor manufacturing processes, power plants, traffic monitoring equipment, and air conditioning monitoring equipment.

The target process 1 of the rainfall inflow prediction system is an urban stormwater drainage process 1 that includes a flow meter 11, a main flow meter 12, main water level gauges 131-13K, ground rain gauges 141-14M, a radar rain gauge 15, pump well water level gauges 161 and 162, an inlet culvert 17, a stormwater pump well 18, a stormwater pump 19, and an inlet gate 110. The radar rain gauge 15 includes radar rain gauges 1511-15QP for each mesh in a QxP mesh.

Various sensors 11-15, 161, and 162 in the urban stormwater drainage process 1 measure quantities that represent the state of the process and quantities related to operation at predetermined intervals (e.g., 30 seconds, 60 seconds). The values measured by the various sensors 11-15, 161, and 162 are collected and stored by a modular time-series data forecasting device.

The inflow culvert 17 is a culvert into which rainwater sent from sewer pipes and pumping stations flows.
The stormwater pump well 18 is a water tank that stores the rainwater that flows in from the inflow conduit 17. The rainwater that flows into the stormwater pump well 18 may have sand and other contaminants removed in advance in a grit basin.
The rainwater pump 19 is a pump that forcibly pumps the rainwater stored in the rainwater pump well 18 into a river or the like.

The inflow gate 110 is installed in the flow path between the inflow culvert 17 and the stormwater pump well 18, and operates to adjust the amount of stormwater flowing from the inflow culvert 17 to the stormwater pump well 18 by opening and closing it. Water level gauges 161 and 162 are installed before and after the inflow gate 110, and the operation of the inflow gate 110 is controlled, for example, according to the values measured by the water level gauges 161 and 162.

FIG. 2 is a diagram illustrating an example of the configuration of a modular time-series data prediction device according to an embodiment.
The modular time-series data prediction device of this embodiment includes a data collection and storage unit 2, a data extraction unit 3, a pairwise prediction model identification unit 4, a prediction model synthesis method definition unit 5, a pairwise output variable prediction unit 6, a synthesis output variable prediction unit 7, a prediction error evaluation unit 8, a pairwise prediction model correction unit 9, and an output prediction result observation unit 10.
The configuration of the modular time-series data prediction device of this embodiment is, for example, a computing device including at least one processor such as a CPU or an MPU, and a memory in which a program executed by the processor is stored, and various functions can be realized by software or a combination of software and hardware.

The data collection and storage unit 2 collects and stores the values of process variables measured by various sensors 11-15, 161, and 162. The data collection and storage unit 2 includes an output variable data selection unit 21 and an input variable data selection unit 22.

The output variable data selection unit 21 selects at least one of the quantities measured by the various sensors 11-15, 161, and 162 as an output variable.
The input variable data selection unit 22 selects a plurality of input variables from the values measured by the various sensors 11-15, 161, and 162. The number of the plurality of input variables may be equal to or less than the number of the plurality of process variables.

The data extraction unit 3 extracts predetermined data from the time-series data stored in the data collection and storage unit 2. The data extraction unit 3 includes an offline prediction model identification data extraction unit 31, an online prediction data extraction unit 32, and an evaluation data extraction unit 33.

The offline prediction model identification data extraction unit 31 extracts data measured by the various process sensors 11-15, 161, and 162 for a specified period (prior to the specified time) from the data collection and storage unit 2 at a specified identification interval or in response to an external request, as data for prediction model identification.

The online prediction data extraction unit 32 extracts prediction data required for online prediction from the measurement values of various process sensors 11-15, 161, and 162 in real time from the data collection and storage unit 2 at a predetermined prediction period.

The evaluation data extraction unit 33 extracts, in real time, the measurement values of the various process sensors 11-15, 161, and 162 required for prediction error evaluation from the data collection and storage unit 2 at a predetermined evaluation interval or in response to an external request.

The pairwise prediction model identification unit 4 uses the measurement value data extracted by the offline prediction model identification data extraction unit 31 to pair the output variables selected by the output variable data selection unit 21 and the input variables selected by the input variable data selection unit 22 for each variable, and identifies a prediction model using the method described below.

The prediction model synthesis method definition unit 5 defines a method for synthesizing output variables (for each output variable if there are multiple output variables) from the output values (pairwise prediction values) of the pairwise prediction model identified by the pairwise prediction model identification unit 4.

The pairwise output variable prediction unit 6 inputs prediction data corresponding to each input variable extracted from the online prediction data extraction unit 32 into the pairwise prediction model identified by the pairwise prediction model identification unit 4 at a predetermined monitoring cycle in the time progression direction, and performs pairwise predictions of the output variable for each input variable (for each output variable if there are multiple output variables).

The composite output variable prediction unit 7 outputs a predicted value of the output variable synthesized according to the method defined in the prediction model synthesis method definition unit 5, based on the pairwise predicted value output by the pairwise output variable prediction unit 6.

The prediction error evaluation unit 8 stores a prediction value that is a composite of the pairwise prediction value output by the pairwise output variable prediction unit 6 and the pairwise prediction value output by the composite output variable prediction unit 7, and evaluates the error between the actual value of the output variable extracted by the evaluation data extraction unit 33 and the stored pairwise and composite prediction values at a predetermined interval or at a predetermined timing.

The pairwise prediction model correction unit 9 determines whether or not adjustment of the pairwise prediction model is necessary based on the prediction error evaluated by the prediction error evaluation unit 8, and performs the necessary pairwise prediction model correction. The pairwise prediction model correction unit 9 can, for example, adjust the parameters of the pairwise prediction model, delete the pairwise prediction model, or integrate the pairwise prediction models. The pairwise prediction model correction unit 9 supplies the corrected pairwise prediction model to the pairwise output variable prediction unit 6.

The output prediction result observation unit 10 is equipped with a display means such as a monitor and a user interface, and is capable of presenting at least the predicted value of the composite output variable by the composite output variable prediction unit 7 to the observer. The output prediction result observation unit 10 may be, for example, a personal computer, a tablet terminal, or a mobile terminal such as a smartphone. If necessary, the output prediction result observation unit 10 may further present to the observer at least one of multiple outputs of the pairwise prediction model by the pairwise output variable prediction unit 6 and the predictive distribution information synthesized from them.

Next, the operation of the modular time-series data prediction device of this embodiment will be described.
For example, in an urban stormwater drainage process 1, various process sensors 11-15, 161, and 162 measure process information at predetermined intervals, and the measured values (values of process variables) are supplied to the data collection and storage unit 2. The data collection and storage unit 2 stores the collected measured values as time-series data in accordance with a predetermined format.

The output variable data selection unit 21 selects variables to be predicted from among the process variables measured by the various process sensors 11-15, 161, and 162, and sets the selected variables as output variables.
For example, in the urban stormwater drainage process 1, when the inflow volume into the normal stormwater pump well 18 is to be predicted, the output variable data selection unit 21 selects the stormwater pump well inflow volume measured by the flow meter (stormwater pump well inflow volume meter) 11 as the output variable to be predicted.

Also, if a flow meter (rainwater pump well inflow meter) 11 is not installed, the inflow volume into the rainwater pump well 18 can be calculated from the measurement values of the rainwater pump well water level gauges 161, 162 and the operational status of the rainwater pump 19 and inflow gate 110, and this can be selected as the prediction target. In this case, the output variable data selection unit 21 calculates and saves the calculated value of the inflow volume, and selects the calculated value of the inflow volume as the output variable.

Furthermore, for example, when predicting multiple items in the urban stormwater drainage process 1, the output variable data selection unit 21 may select multiple items as output variables. Even in the urban stormwater drainage process 1, for example, when stormwater flows into the stormwater pump well 18 from multiple locations, the output variable data selection unit 21 may select the inflow volume at each of the multiple inflow locations as an output variable.

The input variable data selection unit 22 selects as candidate input variables measured values that may affect the output variables selected by the output variable data selection unit 21. If there is no prior information as to whether or not a variable will have an effect, the input variable data selection unit 22 may automatically set all process variables as input variables. The input variable data selection unit 22 may also select variables selected as output variables as input variables. For example, in time series data analysis, so-called autoregression is often considered, in which future values of the variable selected as the output variable are predicted from its own past values. In such cases, the stormwater pump well inflow rate measured by the flow meter (stormwater pump well inflow rate meter) 11 may be selected as both an output variable and an input variable.

Furthermore, when there are multiple output variables selected by the output variable data selection unit 21, the input variable data selection unit 22 may select another output variable as an input variable for a given output variable. This makes it possible to consider the correlation between multiple output variables. On the other hand, when multiple output variables are inflow amounts at different inflow points, there is often no correlation between the output variables due to the physical structure of the process. In this case, the input variable data selection unit 22 can also exclude another output variable as an input variable for a given output variable.

The input variable data selection unit 22 can also exclude from the input variables in advance other variables, such as variables that are clearly unrelated to the output variables based on information such as the physical connection relationships of pipes, or variables for which the causal relationship (the relationship between input (cause) and output (result)) does not physically hold.

For example, it is physically clear that the rainwater pump well water level measured by the rainwater pump well water level gauges 161 and 162 will change due to the influence of the rainwater pump well inflow measured by the flow meter (rainwater pump well inflow meter) 11 selected as the output variable, so the rainwater pump well water level is excluded from the input variables.

The input variable data selection unit 22 may perform a preliminary correlation analysis between variables with unclear physical structural relationships and the output variables selected in advance by the output variable data selection unit 21, and may exclude variables with absolute values of correlation coefficients less than a predetermined value (e.g., 0.5) from the input variables in advance. Furthermore, when performing the correlation analysis, the input variable data selection unit 22 may perform the correlation analysis while shifting the time at which the variable being analyzed is measured, determine the shifted time L at which the correlation is highest, and, based on the sign information of time L, exclude variables with signs that do not satisfy a causal relationship from the input variables in advance. It is likely that the stormwater pump well water level variable mentioned above will also be excluded from the input variables through such an analysis.

In the modular time series data prediction device of this embodiment, the input variable data selection unit 22 selects as input variables the rainwater pump well inflow measured by the flow meter (rainwater pump well inflow meter) 11, excluding only the rainwater pump well water level measured by the rainwater level gauges 161 and 162, the main inflow measured by the main flow meter 12, the main water level measured by the main water level gauges 131-13K, the ground rainfall measured by the ground rain gauges 141-14M, and the radar rainfall at each of the meshes 1511-15QP of the QxP mesh.

In the modular time series data prediction device of this embodiment, for example, K main water level gauges 131-13K are often installed with a relatively large number of main water level gauges, with K being around 10-100, due to recent government trends toward accurately grasping water level information in sewer pipes. Furthermore, weather radar can be installed in many locations because the Ministry of Land, Infrastructure, Transport and Tourism has installed a weather radar called XRAIN nationwide. Depending on the size of the target stormwater drainage district, it is not uncommon for the mesh size of a QxP mesh to be around 100. Therefore, the number of time series data items (number of variables) actually measured is quite large. Furthermore, the number of measured variables typically ranges from hundreds to thousands, not just in urban stormwater drainage process 1, to which the modular time series data prediction device of this embodiment is applied.

In the modular time series data prediction device of this embodiment, the input variables are, for example, the rainwater pump well inflow measured by the flow meter (rainwater pump well inflow meter) 11, excluding only the measurements by the rainwater pump well level gauges 161 and 162; the main inflow measured by the main flow meter 12; K main water level measurements by the main water level gauges 131-13K; M ground rainfall measurements by the ground rain gauges 141-14M; and the radar rainfall measurements for each mesh of the Q×P mesh 1511-15QP. In the following explanation, the number of input variables is assumed to be p (a positive integer).

The offline prediction model identification data extraction unit 31 first acquires time series data for the specified period from the output variable data selection unit 21 and the input variable data selection unit 22 at the requested timing using data for the specified period specified by the person building the prediction model. If the offline prediction model identification data extraction unit 31 periodically builds a prediction model and periodically updates the model, it may periodically extract time series data of input and output variables for the specified period at a specified cycle TL in the direction of time progression. The dataset extracted in this manner is represented as uk, where k = 1, 2, ... p. uk is a column vector, with each row corresponding to a sample at each time in the time series data. Similarly, the time series dataset for stormwater inflow, which is the output variable, is assumed to be y. The offline prediction model identification data extraction unit 31 extracts dataset uk and dataset y at a specified cycle.

Next, the pairwise prediction model identification unit 4 uses uk and y to identify a prediction model for each of the p pairs (u1, y), (u2, y), ... (up, y). In the modular time series data prediction device of this embodiment, the output variable y is only one, the stormwater pump well inflow rate measured by the flow meter (stormwater pump well inflow rate meter) 11. However, if there are multiple (say h) outputs y, the same processing can be performed for each of the h y's. Even if the number of output variables increases to l, the same processing can be repeated h times. Therefore, without loss of generality, we will assume that there is only one output variable.

The most typical prediction model is one that performs identification by assuming a linear regression model that can be expressed by the following equation:
y(t)=a1×uk(tL)+a2×uk(tL-1)+…+an×uk(tL-n+1)+c (5)
Here, a, k = 1, 2, ..., n and c are parameters to be identified. c is a bias parameter that represents the difference between the averages of the input uk and the output y. However, by removing the average values from uk and y in advance and setting the average to 0, c can always be set to 0. Therefore, for simplicity, c will be set to 0 below.

In the fields of digital signal processing and system identification, when uk in equation (5) above is not y itself, it is called a finite impulse response (FIR) model (it can also be called a weighted moving average model for input variables). On the other hand, when uk is equal to y itself, it is called an autoregressive (AR) model. What is important here is that even though it is in the form of linear regression in equation (5), it is not a multi-input, single-output (multi-input, multi-output when multiple outputs are considered simultaneously) linear regression form (such as ARX (Autoregressive eXogenous Input)) in which finite impulse responses from multiple input variables are added together and then autoregression is added, but rather it is still a linear regression model with one input and one output; this is the meaning of "pairwise."

The first advantage of pairwise identification is that, because the concept of "multiple inputs" itself is eliminated, the problem of multicollinearity mentioned earlier (which arises from the relationship between multiple input variables) never arises. This essentially avoids this problem, making it more advantageous than multiple-input, single-output systems in terms of parameter identifiability. To expand on this a bit, in the field of system identification, the condition for parameter identifiability (i.e., parameters can be uniquely determined) known as the persistent excitation condition (PE condition) is that the number of identifiable parameters changes depending on the amount (number) of frequency components contained in the input signal uk. It is known that the greater the number of parameters ak to be identified (estimated), the greater the number of frequency components the input signal uk must contain.

On the other hand, in the modular time series data prediction device of this embodiment, the values of the input variables cannot be controlled (adjusted) by the builder of the prediction model, but are given observation (measurement) information. Therefore, the number of frequency components included is given as an external condition. The fewer the number of parameters, the better the identifiability, and conversely, the more parameters there are, the greater the likelihood that identifiability will deteriorate. More precisely, when the number of parameters is small, the likelihood of identifiability deteriorating is 0%, and identifiability with a small number of parameters is higher, or at worst, the same as with a large number of parameters. Therefore, while the one-input, one-output model in equation (5) above contains n parameters, a p-input, one-output prediction model will contain p x n parameters. Therefore, it can be seen that the greater the number of inputs p, the greater the likelihood that it will be difficult to uniquely identify the parameters.

Furthermore, if there is a dependency (correlation) between the input variables mentioned above, the problem of multicollinearity arises, and identifiability is lost regardless of the number of parameters to be identified. Because the input information is given, it is, in principle, impossible to uniquely determine parameter values no matter what method is used. In contrast, the problem of multicollinearity does not, in principle, arise in the one-input, one-output model of equation (5). The only time parameters cannot be identified in the form of equation (5) is when the input uk does not contain sufficient frequency information. In this case, by adjusting the order n and the number of parameters to be identified, it is possible to ensure that parameter ak is identifiable (uniquely identifiable).

For example, if the value of uk is constant, then uk(t-L) = uk(t-L-1) = ... = uk(t-L-n+1). In such cases, parameter identifiability can be maintained by setting n = 1 and identifying only one parameter, a1. Generally, the upper limit nmax of the number of identifiable parameters n can be determined by examining the above continuous excitation conditions using a given uk. Therefore, by using various criteria and standards, such as AIC (Akaike Information Criterion), BIC (Bayesian Information Criterion), and MDL criterion (Minimum Decoding Criterion), or by performing cross-validation within the range of n≦nmax, an appropriate value for n can be determined, thereby consistently maintaining the identifiability of equation (5).

The second advantage of pairwise identification is that even if there is only one input variable excluding the output variable, this method separates and identifies two models: an autoregressive AR model that uses the output variable as input, and an FIR model that uses the input variable (an input variable that is not an output variable) as input. In such cases, conventional non-pairwise methods would simultaneously identify the autoregressive component (AR component) and weighted moving average component (FIR component/MA component) using an ARX model or the like. However, in many cases, the correlation between the output variable itself (autocorrelation) has a stronger effect than the relationship (correlation) between the output variable and the input variable, and the autoregressive component often has an extremely strong effect on the prediction. In such cases, if the value of the input variable changes suddenly, the system will be unable to track the change, and predictions are likely to be delayed.

In the case of the modular time series data forecasting device of this embodiment, for example, if the output variable is pump well inflow and the input variable is ground rainfall at one location, the influence of the autoregressive component on the inflow will be stronger than the influence of the FIR component on the ground rainfall. Even if the ground rainfall value changes suddenly due to a change in rainfall, the change in the predicted value will be small until a change in the inflow appears, and the predicted value will likely lag behind the actual measured value. In contrast, with a pairwise method that separately identifies an autoregressive model (AR model) for the inflow and a finite impulse response model (FIR model) for the ground rainfall, a user such as a supervisor can adjust the influence of the AR component and the FIR component on the forecast.

The third advantage of pairwise identification, which is closely related to the second advantage, is that the meaning of the model in equation (1) for each input variable is clear. Here, as before, consider the simple ARX model in equation (2) below, which has one input variable excluding the output variable.
y(t)=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny)
+b1×u(tL)+ b2×u(tL-1)+…+bnu×u(tL-nu+1) (6)

In this case, if the above equation (6) is identified as an ARX model, the meanings of the following terms related to y (defined as Y) and u (defined as U) are unclear, and it is difficult to explain the individual meanings of Y and U.

Y:=a1×y(t-1)+ a2×y(t-2)+…+any×y(t-ny) (7)
U:=b1×u(tL)+ b2×u(tL-1)+…+bnu×u(tL-nu+1) (8)
On the other hand, since equations (7) and (8) each have the form of equation (5), when pairwise identification is performed, Y and U each mean a predicted value of y, and therefore Y and U can be interpreted as a predicted value of output y for each separate explanatory variable (input variable).

By performing pairwise identification in this way, it is easy to construct a predictive model in the form of equation (6) constructively from equations (7) and (8). For example, by averaging U in equation (7) and Y in equation (8) (= (U + Y) ÷ 2), a predictive model in the form of equation (6) can be obtained, which clearly again represents the predicted value of y.

As mentioned above, even if the form (structure) of a predictive model is the same, changing the identification procedure will change its meaning. The advantage of pairwise identification is that it allows for the construction of a model that is easy to interpret.
These are the advantages of identifying a prediction model pairwise with one input and one output, rather than with multiple inputs and one output. In addition, by identifying the pairwise prediction model in the form of an FIR or AR regression equation as in equation (5), various analysis methods well known in the fields of digital signal processing and system identification, such as frequency analysis and stability analysis, can be directly applied to pairwise analysis. Therefore, if you have knowledge in that field, you can perform more detailed analysis and adjustments in pairwise analysis.

As a pairwise prediction model, a prediction model can be constructed by further decomposing an AR model/FIR model such as equation (5) into a simple regression model in which each variable taking time delays into account is considered as a single explanatory variable.
y(t)=a×uk(tL)+c (9)
When simplified as in equation (9) above, the parameters included in equation (9) are only three: delay time L, the magnification (gain, proportionality coefficient) a of the output variable relative to the input variable uk, and the difference between the average values of the input variable and the output variable (output value when the input is 0, bias) c, and the meanings of all parameters become clear.

In equation (5), the interpretation of delay time L and bias c is the same as in equation (9), but the individual meanings of the regression coefficients (parameters) a, k = 1, 2, ..., n are unclear. While equation (5) itself can be interpreted using knowledge of system identification and digital signal processing, intuitive interpretation of the individual a is difficult. In contrast, in the simple regression equation of equation (9), a is easy to interpret; it simply indicates the magnification of the output variable relative to the input variable when the mean value is set to 0. This makes it easy to interpret the meaning of all parameters. Therefore, using equation (9) makes it easy for predictive model builders to intuitively and manually adjust parameters when the parameter identification results from applying an algorithm deviate from intuitive perception (such as can occur when anomalous data is mixed in), and it also makes it extremely easy to explain the contents of a predictive model to others.

In this embodiment of the modular time series data forecasting device, for example, if the output variable is the pump well inflow and the input variable is ground rainfall at one location, then c in equation (9) represents the inflow when ground rainfall is zero. In the case of a combined sewer system in which sewage and rainwater pass through the same pipe, this can be interpreted as equivalent to the sewage volume. Furthermore, a in equation (9) represents the magnification of the rainwater inflow volume relative to the rainfall measured by the ground rain gauge when the sewage volume is subtracted from the inflow volume. If the correlation between ground rainfall and inflow volume is sufficiently high, even this simple magnification is meaningful. Furthermore, the delay time L in equation (9) represents the delay time for rainwater flow from the ground rainfall observation point to the stormwater pump well, and can be interpreted as the sum of the outflow delay before it flows into the sewer pipe and the flow time through the sewer pipe to reach the stormwater pump well.

In this way, the advantage of using a simple model such as equation (9) is that it is extremely easy to interpret and explain. In fact, the simplest model using only equation (9) is also adopted in Patent Document 1 as a simplified model for inflow volume prediction.

A specific example will be described in which the above-described intuitive interpretation makes it easier to manually adjust parameters.
3 and 4 are scatter diagrams that schematically show an example of the relationship between main water level data and inflow data.

Figures 3 and 4 show an example of a scatter plot in which the horizontal axis represents main water level data at a location corresponding to one of the main water level gauges 131 to 13K, and the vertical axis represents inflow data to be predicted corresponding to inflow meter 11.
The straight line shown in Figure 3 is the gain and bias of simple regression calculated by the least squares method, and in this case the gain a = 9.6024 and the bias c = 23.3142. This is the optimal regression line in the sense that it minimizes the mean square error (MSE or RMSE) of the predicted value.

In Figure 3, for example, when the main water level is -1.5, the regression line predicts that the inflow value will be approximately 9. However, the scatter plot in Figure 3 shows that when the main water level is -1.5, the inflow value is overwhelmingly greater than 9.

As mentioned above, the reason for the discrepancy between the predicted value based on the regression line and the measured value is that there is an overwhelming amount of data for main water levels around -2.5 to -2, so the coefficients are estimated to better fit the data in this range. While this regression line is correct in theory, it is often not necessarily desirable for actual inflow predictions. This is because the predicted inflow values are used to control and assist the operation of stormwater drainage pumps, and in order to avoid the risk of flooding, there is an implicit desire to predict cases where inflow volumes are high as accurately as possible, and prediction accuracy when inflow volumes are low is often not as important.

However, incorporating such implicit requirements into an algorithm requires measures such as weighting the prediction errors in areas with high inflows, which generally requires a considerable amount of effort, including adjusting those weights.

Figure 4 shows an example of regression lines L1-L11 obtained by applying the least squares method and 10 other regression methods to the same data as in Figure 3. Of these, methods 1 through 10 are various advanced regression methods built into a software package called Python, but none of the regression lines L1-L10 obtained by any of these methods were able to produce appropriate predictions when the inflow volume was large.

The 11th method extracts only data from the main water level range of -2.2 to -1.5, and applies the ordinary least squares method without using any other data. The inventors of this application, understanding that a regression line can be obtained to fit data for main water levels around -2.5 to -2, intentionally discarded data in this area, and intentionally adjusted the method to fit data when inflow volumes are high.

As such, when attempting to perform adjustments that rely on algorithms, it may be necessary to take various measures, such as limiting the range of data or thinning out the data, to obtain the desired parameters. However, if the straight line L11 obtained by the 11th method is based on the extremely simple principle that a straight line can be uniquely determined by determining two points, it takes almost no effort for a human to visually draw a line like the straight line L11. As described above, by using simple regression, which is extremely easy to understand, and parameters that can be clearly visualized on the scatter diagrams shown in Figures 3 and 4, it can be seen that manual adjustments can be made extremely easily.
In other words, when a simple regression model is used, it can be interpreted using easily understandable parameters, namely, gain and bias, making adjustments easy.

Furthermore, simple regression allows for interpretation from a different perspective. That is, it allows for interpretation using three parameters, the mean, standard deviation, and correlation, which are basic statistical quantities in the field of statistics.

It is widely known that the regression coefficient (proportional coefficient, slope, magnification, gain) a of a simple regression model has the following relationship with the correlation coefficient r:
a = r × (σy/σu) (10)
where r, σy, and σu are the correlation coefficient between the input variable and the output variable, the standard deviation of the output variable, and the standard deviation of the input variable, respectively. It is also widely known that the regression coefficient a and bias c have the following relationship when identified by the least squares method so as to minimize the prediction square error of the simple regression model.
c = μy - a × μu (11)
Here, μy and μu represent the mean of the output variables and the mean of the input variables, respectively.

Therefore, in equation (9), if the input variables and output variables are normalized using the mean and standard deviation, respectively, the values of the regression coefficient a and the correlation coefficient r will match, and c=0.
This relationship means that explaining the model using the gain and bias of a simple regression model is equivalent to explaining the model using only the mean and standard deviation of each variable (input variable and output variable) and the correlation coefficient between the input and output. Therefore, if you normalize the data for each variable in advance using appropriate estimates of the mean and standard deviation for each variable, the input-output relationship of simple regression can be explained entirely using only the correlation coefficient, as shown in the following equation.
y´(t)=a×uk´(tL)=r×uk´(tL) (12)
where y' = (y - μy)/σy and uk' = (uk - μuk)/σuk are the normalized output and input variables.

When interpreted in this way, for example, if there is no correlation between input and output, i.e., if the input-output relationship is such that r = 0, then the input corresponding to that output has no predictive ability whatsoever, as can be clearly seen from the relationship in equation (11). Furthermore, real-world data is often contaminated with various noises and outliers. In such cases, a robust estimation method that robustly estimates the three parameters (mean, standard deviation, and correlation coefficient) against outliers can be used to construct a predictive model that is robust against outliers. Of course, various robust estimation methods have been developed for estimating the parameters of more complex models, but the results of robust estimation are generally not easy to interpret. However, robust estimation of the three parameters (mean, standard deviation, and correlation coefficient) makes it easy to interpret and judge the quality of the results of any robust estimation applied. This also facilitates interpretation and adjustment when building a predictive model for data contaminated by noise or outliers.

As such, equations (9) and (10) through (12) are composed only of easily understandable parameters. However, compared to equation (5), they have an extremely simple structure, and in actual applications, their expressive power may be insufficient, resulting in insufficient prediction accuracy. Therefore, we consider synthesizing equation (5) from equation (9) (or equations (10) through (12)). To obtain equation (5) using the model of equation (9), we simply make the delay time L in equation (9) variable and add them together. That is, first, we match the minimum value of the delay time L in equation (9) (hereafter referred to as L' to distinguish it from L in equation (5)) to L in equation (5), setting L' = L. Then, for delay times greater than L, add n models of equation (9) where L' = L+1, L' = L+2, ..., L' = L+n-1 to obtain a model in the form of equation (5).

That is, by adding together n simple regression models of the following equation (13), a model in the form of an FIR or AR model of equation (14) which is formally equivalent to equation (5) can be obtained.
y(t)=a1×uk(tL)+c (13_1)
y(t)=a2×uk(tL-1)+c (13_2)
y(t)=an×uk(tL-n+1)+c (13_n)
y(t)=1/n×(a1×uk(tL)+a2×uk(tL-1)+…+an×uk(tL-n+1))+c (14)

By creating n simple regression models for n different delay times and then creating an average model of these, a model in the form of an FIR or AR model, formally the same as equation (1), can be obtained. In this case, a simple averaging process can be performed as above, but a weighted average can also be performed as shown below, using the relationship between the gain, bias, and correlation coefficient of the simple regression described above.

First, input and output variables are normalized in advance to express equation (9) in the form of equation (12). Then, the regression coefficient coincides with the correlation coefficient, so an equation in the form of an FIR/AR model, equation (15) similar to equation (14), can be obtained as a weighted average weighted by (the absolute value of) this regression coefficient (hereinafter, no distinction will be made between the symbols for normalized and non-normalized variables, and u(t) and y(t) will represent normalized variables depending on the context).
y(t)=(a1×uk(tL)+…+an×uk(tL-n+1))/(|a1|+|a2|+…+|an|)) (15)

This makes it possible to synthesize pairwise prediction models that are weighted for predictions based on lag times where correlation is stronger. Even if the relationship between the regression coefficient and the correlation coefficient cannot be explicitly determined as in equation (12), it is possible to perform weighted averaging using the absolute value of the correlation coefficient or the absolute value of the regression coefficient. However, the existence of the relationship in equation (12) makes the significance of weighting clearer and more persuasive. Furthermore, in equation (15), weighting can also be performed using the square of the correlation coefficient rather than the absolute value of the correlation coefficient.

As pairwise prediction models, both of the above synthesis methods can only express linear relationships between input and output, but there are often nonlinear relationships between input and output. In such cases, instead of equation (5), an appropriate nonlinear function φ(·) can be used to perform a nonlinear transformation into a nonlinear regression form such as equation (16) below.

y(t)=a1×φ(uk(tL))+…+an×φ(uk(tL-n+1))+c (16)
This is a nonlinear regression technique that is widely known in the field of machine learning, and various linear regression methods can be directly applied as long as linearity with respect to the parameters a, k = 1, 2, ..., n is maintained. Also, as is well known in the field of machine learning, instead of directly specifying the nonlinear function φ(·), it is possible to directly specify a kernel function (similarity function) that can be converted to be equivalent to specifying φ(·) for a matrix called a design matrix that is required when calculating equation (16). However, in this case, it should be noted that intuitive interpretability may be slightly reduced.

By using nonlinear regression such as equation (16), it is possible to handle cases where there is a nonlinear relationship between a certain input variable and an output variable. In identifying pairwise prediction models, a regression model with the output variable is constructed for each input variable. Therefore, if the relationship between one input variable and the output variable is nonlinear but the relationship between another input variable and the output variable is linear, it is sufficient to simply use equations (5) and (16) for each input variable, which is generally expected to improve interpretability. The above series of actions are the actions of the pairwise prediction model identification unit 4 in this embodiment.

Next, the prediction model synthesis method definition unit 5 defines a synthesis method for the prediction models defined by the pairwise prediction model identification unit 4 .
The simplest method for combining prediction models is to combine prediction models by averaging the p prediction outputs (hereinafter, the prediction outputs for the input variables u1, u2, ..., up) of the pairwise prediction model identification unit 4, as follows:
y(t)=mean(y1(t), y2(t),..., yp(t))=1/p×(y1(t)+y2(t)+...+yp(t)) (17)
The definition of equation (17) is the most basic definition method, but by improving this definition, it becomes possible to improve the reliability of the synthesized prediction output and to perform prediction with an intentional tendency (bias) in the prediction. This will be explained in order below.

First, when performing a simple averaging process (sample mean) such as in equation (17), the predicted values of each pairwise prediction model are treated equally. This means that predicted values from pairwise prediction models with good prediction accuracy may be mixed with predicted values from pairwise prediction models with poor prediction accuracy, potentially degrading the overall prediction accuracy of the combined model. Furthermore, in actual operation, it is not uncommon for measurement data for one input variable to be unreliable while measurement data for another input variable is highly reliable. However, averaging the predicted values of pairwise prediction models for both highly reliable and unreliable input variables can degrade the prediction accuracy of the combined prediction output. In extreme cases, such as when a sensor failure or malfunction for a certain input variable causes the time series data for that input variable to contain a large number of outliers (abnormal values), not only does prediction accuracy deteriorate, but the outliers can also render the combined prediction meaningless (even causing it to fail).

To deal with such a situation, the prediction model synthesis method definition unit 5 can take the following multiple approaches.
The first is an approach that places emphasis on ensuring that the synthesized predicted value itself does not break down, mainly for the latter outliers, and is a method of adopting robust estimation instead of the sample averaging process of equation (17).

The "sample mean" process used in equation (17) is one method of estimating a representative value from the p predicted values y1, y2, ..., yp, and in the field of statistics, this type of representative value estimation is called location parameter estimation. In other words, the "sample mean" is one method of estimating a location parameter, and statistically, the sample mean is known to have good properties in terms of estimation efficiency, but it is also known to be the most vulnerable (non-robust) in terms of robustness against outliers. Therefore, in order to improve robustness against outliers, it is better to replace equation (17) with a robust location parameter estimation.

In the field of robust estimation, several indices are known for evaluating robustness against outliers, but the one that is easiest to understand intuitively is the breakdown point. This is an indicator of what percentage of outliers is acceptable to be mixed into the data used in statistical estimation. If X% of the data used is replaced with ∞ to assume hypothetical outliers, the value of X at the boundary (just before ∞) where the estimator (the estimated value of the location parameter, when estimating a location parameter such as the mean) reaches ∞ (i.e. breaks down) or not is called the breakdown point.

For example, when estimating the location parameter known as the "sample mean," simply replacing one point of data with ∞ also makes the sample mean value ∞, so the breakdown point of the "sample mean" is 0%. In robust statistics, it is known that the maximum breakdown point is 50%, and the "median" is known as the location parameter estimator with the largest breakdown point. Therefore, if the primary goal is to ensure that predicted values are not affected by outliers, the "mean" in equation (17) can be replaced with the "median" to form the definition in the prediction model synthesis definition unit 5.

On the other hand, while the "median" is robust, it is known to have poor estimation efficiency. It is intuitively obvious that, since it only uses one predicted value out of p variables, it is not the best method when high-accuracy predictions are important. For this reason, various robust estimation methods are known that achieve both improved estimation efficiency and robustness. The simplest method is a process called the "trimmed mean," which removes α% of the top and bottom values of p data and then takes the average. As is intuitively clear, increasing α improves robustness, with median estimation as the limit. Reducing α toward 0 clearly decreases robustness, with the limit being the normal mean (sample mean). Therefore, this trimmed mean can also be used as the definition in the prediction model synthesis definition unit 5. In this case, robustness can be adjusted by adjusting α, but since it may be difficult to adjust it appropriately, another robust estimation method can be used.

A typical example of such a method is to perform location parameter estimation known as the Hodges-Lehmann estimator (HL estimation), which calculates the average of two yi and yj (i ≠ j) for all combinations (p(p-1)/2) of p values y1, y2, ..., yp, and then uses the median. This method is known to be highly robust with a breakdown point of approximately 30%, as well as having good estimation efficiency, and this HL estimation can also be defined in the prediction model synthesis method definition unit 5.

Another robust estimation method known as MCD, which is generally used for multivariate data, can also be applied. MCD extracts a specified percentage (usually 50% to 75%) of p data items and performs estimation using the data with the smallest variance among all combinations. This method can also be used to estimate the mean. MCD is also known as a highly robust and efficient estimation method. On the other hand, HL estimation and MCD estimation require calculations for all combinations of the number of data items extracted from p data items (2 for HL estimation, approximately p/2 to 3p/4 for MCD). Therefore, as the number of input variables p in the pairwise model increases, processing time increases dramatically. Therefore, even if they are used as the definition for the prediction model synthesis method definition unit 5, real-time performance may not be ensured when actually making predictions.

To address such a situation, instead of processing "all combinations," a bootstrap method may be used to repeatedly perform processing as many times as possible in real time. That is, a predetermined k number of data items are randomly and repeatedly extracted from p data items y1, y2, ..., yp, and the average value is calculated a predetermined number of times (N times). From the N average values calculated repeatedly, for example, the median of the N average values is used, as in HL estimation. By defining such average value estimation using the bootstrap method in the prediction model synthesis method definition unit 5, it is believed that it is possible to synthesize prediction values with high estimation efficiency and robustness while maintaining real-time performance.
In addition to the above, the prediction model synthesis method definition unit 5 can also define a method in which a robust estimation such as M estimation, which performs estimation by weighting outliers, is applied to the location parameter estimation.

The second approach aims to perform a composite prediction by weighting predicted values with good prediction accuracy, whereas the previous approach aimed primarily at estimating a representative value of predicted values robustly against outliers (outlier predicted values). The basic idea is to perform a weighted average instead of the simple sample average in equation (17), and the composite prediction output is defined by the following equation:
y(t)=w1×y1(t)+w2×y2(t)+…+wp×yp(t) (18)
Here, wk, k = 1, 2, ..., p are weights that satisfy the constraint that w1 + w2 + ... + wp = 1. In the definition of the prediction model synthesis method definition unit 5, it is necessary to define a method for setting these weights.

For this purpose, the pairwise prediction model identification unit 4 can use the (multiple) correlation coefficient between the predicted output when the pairwise prediction model is identified and the actual output, or the coefficient of determination, which is its square. This makes it possible to synthesize predicted outputs using a weighted average that places emphasis on the predicted output of models with high prediction ability and high prediction accuracy, at least for the data used when identifying the pairwise prediction model.

Furthermore, if a pairwise prediction model is constructed, for example, as in equation (15) above, the correlation coefficients between the explanatory variables and output variables for each delay time have already been calculated, so for each integrated pairwise prediction model (corresponding to equation (14)), the average value rk, k = 1, 2, 3, ..., p of the absolute values of each correlation coefficient can be calculated, and the weight wk can be defined as wk = rk/(r1 + r2 + ... + rp) according to rk.

In this way, the operation of determining the weights of the weighted average values according to the predictive capabilities of the pairwise prediction models and defining the prediction model synthesis method using equation (18) is an example of the operation of the prediction model synthesis method definition unit 5 in this embodiment.

The third method, like the second method, defines the combined prediction output using the weighted average of each pairwise prediction model in equation (17), but determines the weights by identifying (estimating) them using data. In this case, by using the data used when identifying the pairwise prediction model and the identified parameter values, time series data of the predicted values y1, y2, ..., yp of each pairwise prediction model for the identified data can be obtained.

If we consider equation (17) with the predicted values y1, y2, ..., yp as inputs and the actual output data used for identification as the output y, then equation (17) itself takes the form of a multiple regression model with weights wk, k = 1, 2, ..., p as regression coefficients. Therefore, using this predicted time series data and the output data used for identification, the weights wk, k = 1, 2, ..., p can be identified using multiple regression. However, because the weights must each satisfy the inequality constraint that they must be positive and the equality constraint that w1 + w2 + ... + wp = 1, the least squares method used in ordinary multiple regression cannot be applied. However, it is possible to estimate the weights using a linear regression estimation method that can take constraints into account, such as mixed linear estimation.

This method can also be applied to deriving equation (15) instead of equation (14) from equation (13), so a predictive model can be synthesized through three stages of learning (estimation): simple regression model ⇒ FIR/AR model ⇒ multi-input transfer function model.

In this embodiment, the method of defining the prediction model synthesis method using equation (17) in which the weights wk, k = 1, 2, ..., p are estimated as described above is another example of the operation of the prediction model synthesis method definition unit 5.

The fourth approach, which differs from the first three, is used when you want to impart a bias to the synthesized prediction, and is a method of synthesizing prediction outputs when the goal is to intentionally over- or under-predict.

Situations where the above prediction is required are often encountered in real-world problems. For example, in the case of the rainwater inflow prediction performed in this embodiment, the prediction information is often used as support information for determining the timing of starting and stopping a rainwater drainage pump. In such cases, predicting an inflow amount that is lower than the actual rainwater inflow amount can be fatal from the perspective of risk avoidance.
In such cases, it goes without saying that it is preferable to accurately predict the amount of inflow, but if this is practically difficult, it is safer to predict a slightly larger amount.

Furthermore, for example, when using inflow volume predictions to control a stormwater drainage pump, it is preferable to make an over-prediction to determine the timing to start the pump, but it is preferable to make an under-prediction to determine the timing to stop the pump, so there are cases where it is necessary to make an under-prediction.

Also, for example, when predicting demand for water or electricity, if the prediction is too low, it could cause disruptions to the supply of water or electricity, so it may be safer to predict slightly higher.

Furthermore, by making it possible to make predictions on the over or under side in this way, it becomes possible not only to simply make over or under predictions, but also to make fine adjustments when actually operating the prediction system and performing long-term evaluations, by readjusting the predictions to the over side if they tend to be underestimated, or conversely, by readjusting the predictions to the under side if they tend to be overestimated.

Based on this motivation, the method for synthesizing predicted outputs is described below. To achieve this, in addition to defining the estimation method for the location parameter described above for each pairwise prediction model output predicted value y1, y2, ..., yp, we also simultaneously define estimation methods for the scale parameters corresponding to the standard deviation and variance. For example, if the location parameter is the mean, we define the standard deviation as the scale parameter. Furthermore, if the median (center value) of the robust estimation method is used as the location parameter, we define MAD (Median Absolute Deviation) as the scale parameter. Furthermore, if HL estimation, MAD estimation, bootstrap estimation, M estimation, etc. are used, we define the estimation method for the scale parameter for each location parameter. In this way, we predefine the estimation methods for the location parameters and scale parameters to be applied to y1, y2, ..., yp. Below, we denote these as μ and σ, respectively (μ and σ are usually used to mean the mean and standard deviation, but here we will use them to refer to the location parameter and scale parameter defined here, not necessarily the mean and standard deviation).

Next, for the set of predicted outputs y1, y2, ..., yp of each pairwise prediction model, the value and sign of k are specified so that the predicted output is a combination of μ±kσ values for that set. For example, specifying μ+σ with k=1 means that the predicted output is a combination of values that are one time larger than the scale parameter and one time larger than the location parameter in y1, y2, ..., yp. Conversely, specifying k=-1 means that the predicted output is a combination of values that are one time smaller than the scale parameter and one time larger than the location parameter in y1, y2, ..., yp. If the mean and standard deviation are used as the location and scale parameters, and assuming that y1, y2, ..., yp follow a normal distribution, μ+σ corresponds to extracting approximately 68% of the value (μ is the 50% value), and μ-σ corresponds to extracting approximately 32% of the value.

Incidentally, even without using the mean or standard deviation, you can extract the K% quantile corresponding to K% of y1, y2, ..., yp and use it as the composite output. If the location parameter and scale parameter are defined appropriately, this can be written in the form μ±kσ, so it is also possible to specify it by directly extracting the quantile set in this way. For example, if you want to make an overprediction, you can define a composite forecast output that uses a predicted value corresponding to the 90% quantile to make an overprediction, assuming that you will exclude some outliers in the prediction, or you can define a composite forecast output that uses a predicted value corresponding to the 10% quantile to make an underprediction.

Four methods have been described as definition methods for the prediction model synthesis method definition unit 5, but the prediction model synthesis method definition unit 5 can define a prediction model synthesis method using any of these methods, or a combination of these. The above series of operations is executed offline, either periodically or irregularly, using past data.

Next, prediction is performed online using the definition equations identified and defined offline.
Online prediction is performed at a predetermined period TH (<<TL (period when offline identification is performed periodically)) in the direction of time progression.

First, input variable data is extracted from the online prediction data extraction unit 32 at a period TH. In this embodiment, the input variables include the rainwater pump well inflow rate measured by the flow meter (rainwater pump well inflow rate meter) 11, the main inflow rate measured by the main flow meter 12, K main water level gauges measured by the main water level gauges 131-13K, M ground rainfall rates measured by the ground rain gauges 141-14M, and radar rainfall rates in each of the QxP meshes 1511-15QP.

Next, the pairwise output variable prediction unit 6 calculates p predicted outputs using the pairwise prediction models of each input variable and output variable identified by the pairwise prediction model identification unit 4. This can be calculated immediately by inputting the online input variable data extracted by the online prediction data extraction unit 32 into the pairwise prediction model whose parameter values have been identified in the form of equations (5), (14), (15), etc. This is the function of the pairwise output variable prediction unit 6.

Next, the composite output variable prediction unit 7 uses the definition equations of the synthesis method, which correspond to equations (17) and (18) defined in the prediction model synthesis method definition unit 5, to input the pairwise predicted outputs for each input variable output from the pairwise output variable prediction unit 6, and calculates a synthesized predicted output.

Next, the prediction error evaluation unit 8 obtains the predicted output results calculated by the composite output variable prediction unit 7 and the pairwise predicted output results from the pairwise output variable prediction unit 6, and saves them as time-series data at a predetermined cycle TH in the direction of time progression. The saved prediction results are retained until the time to perform error evaluation. The timing for evaluating the prediction error may be externally instructed by a user of the prediction model, or may be performed at a predetermined cycle TM in the direction of time progression (e.g., TH<TM<TL). In either case, at the time to evaluate the prediction error, the evaluation data extraction unit 33 extracts the output variable measured at the time corresponding to the time-series data of the predicted output saved in the prediction error evaluation unit 8, i.e., in this embodiment, the time-series data of stormwater inflow.

The forecast error evaluation unit 8 then evaluates the error between the time series data of the extracted output (stormwater inflow) and the predicted output of the p pairwise forecast models, as well as the error of the combined forecast output. In this case, the error evaluation can be performed using MSE (mean square error = L2 error), which is typically performed in time series analysis, or its square root, RMSE. However, appropriate error evaluation criteria can be set depending on the purpose. For example, if the prediction accuracy of peak stormwater inflow is important, the error evaluation criteria can be the difference in peak inflow or the difference at the maximum error (L∞ error = ∞ norm error). Evaluation criteria such as the Nash-Sutcliffe coefficient can also be used. Furthermore, if the goal is to improve delays when inflow increases suddenly, the phase difference between the predicted output and actual output can be calculated, and the reduction of the phase difference can be used as the error evaluation criterion. In either case, the forecast error evaluation unit 8 evaluates the error between the time series data and the predicted output of the p pairwise forecast models, as well as the error of the combined forecast output, at a predetermined timing or period TM based on the set error evaluation criteria.

Next, the pairwise prediction model correction unit 9 adjusts the parameters of the pairwise prediction model, deletes specific pairwise models, or integrates multiple pairwise models based on the prediction error evaluated by the prediction error evaluation unit 8.

An example of a more specific method of modular prediction to which the above-described adjustment function, deletion function, and integration function are added will be described below.
For the first parameter adjustment, the pairwise prediction model correction unit 9 extracts input/output pairs that are judged to have significantly poor prediction accuracy (prediction errors larger than a predetermined allowable error) based on the prediction errors evaluated by the prediction error evaluation unit 8, and re-identifies the pairwise prediction models of the extracted input/output pairs. The poor prediction accuracy can be determined by the following methods: (1) directly specifying a threshold (tolerance) for the prediction error and determining a pairwise prediction model (or input/output pair) that exceeds the threshold as having poor prediction accuracy; (2) comparing the error of the combined prediction output with each prediction error of the pairwise prediction model and determining a pairwise prediction model with a prediction error that is significantly larger (exceeding a preset tolerance) than the combined prediction error as having poor prediction accuracy; or (3) calculating the location parameter μe and the scale parameter σe for the prediction error of the pairwise prediction model (which may include the combined prediction error) and determining a pairwise prediction model with a prediction error outside μe±kσe (where k is a parameter that can be set and is usually set to about 2 to 3) as having poor prediction accuracy. If it is determined that parameter readjustment (re-identification) is necessary based on the determinations (1) to (3), the pairwise prediction model correction unit 9 re-identifies the parameters using the most recent data from a predetermined period.

In this embodiment, for example, if it is determined that the accuracy of the pairwise forecast model between any of the M pieces of ground rainfall data and rainwater inflow has deteriorated, the pairwise forecast model correction unit 9 re-identifies only the parameters of the pairwise forecast model between ground rainfall and rainwater inflow. The ability to perform readjustment through partial identification in this way is a major advantage of a modular configuration (configurative approach).

The second specific pairwise model can be deleted using a procedure similar to that used for adjusting the first parameter. The basic procedure involves the pairwise prediction model correction unit 9 evaluating the error using a method similar to that described for adjusting the first parameter, extracting pairwise prediction models with significantly larger errors and degraded accuracy, and deleting the extracted pairwise prediction models. The error evaluation method is the same as the first method, but the pairwise prediction model correction unit 9 extracts pairwise prediction models to be deleted by (1) setting a higher threshold for judging the degree of degradation than the first one, and deleting models whose prediction accuracy is clearly worse than after the first adjustment, or (2) first readjusting the first parameter, and deleting models that still do not meet the criteria used to determine accuracy degradation even after parameter readjustment.

In this embodiment, for example, if a main water level meter becomes submerged, it may become unable to accurately measure water levels, resulting in the inclusion of a large number of outliers (abnormal values). In this case, the pairwise prediction model based on the submerged water level meter will be separated. In this way, a major advantage of a modular configuration (configurative approach) is that prediction accuracy can be improved by partially isolating parts that are negatively affecting prediction accuracy.

The third integration of pairwise prediction models is performed as follows: First, the pairwise prediction model correction unit 9 compares and evaluates the prediction errors of the pairwise prediction models, and extracts and groups pairwise prediction models with similar prediction accuracy, where the absolute value of the difference between their prediction errors is below a predetermined threshold, as similar prediction models.

Next, the pairwise prediction model correction unit 9 compares the parameter values of the pairwise prediction models grouped in this way within the group. For example, the coefficient parameters a k, k = 1, 2, ..., n are vectorized to A = [a1, a2, ..., an], and the similarity between the coefficient parameter vectors of the prediction models being compared is calculated, for example, using the inner product, and pairwise prediction models that are determined to be similar according to a predetermined threshold value are selected as candidates for integration.

In the modular constructive approach of this embodiment, parameter identifiability (unique decidability) is basically guaranteed, so if the prediction results and parameter values are similar, the input variables should have approximately the same values. However, just to be sure, the pairwise prediction model correction unit 9 calculates the correlation coefficient between the input variables of the pairwise prediction models extracted as candidates, and if it is equal to or greater than a predetermined threshold (e.g., a correlation coefficient of 0.95), it determines that the pairwise prediction models should ultimately be merged. The pairwise prediction model correction unit 9 then defines a variable that is the average of the input variables of the multiple pairwise prediction models determined to ultimately be merged as a new, combined input variable.

The pairwise prediction model correction unit 9 then determines the coefficients of the pairwise prediction model for the combined input variables using either (1) adopting the average value of the previously compared coefficient vectors, or (2) re-identifying and estimating the newly combined averaged input and output variables.

In the modular time series data forecasting device of this embodiment, the above-mentioned integration is thought to occur, for example, in the following case. In a case where the input variables include radar rainfall data for each of QxP meshes 1511 to 15QP, the size of the radar rainfall mesh is approximately 250m x 250m for an X-band weather radar called XRAIN, for example, so the rainfall data in adjacent meshes often have almost identical values. In such cases, the rainfall data in adjacent meshes is likely to be integrated, and the QxP input variables are likely to be integrated into a small number of independent composite input variables. When this integration occurs, the radar rainfall mesh data is integrated according to the similarity of its impact on rainwater inflow, which not only reduces the number of input variables but is also expected to improve interpretability.

In addition, in cases like this, the problem of multicollinearity mentioned earlier occurs. To address this issue, if input variables are automatically selected using a method such as Lasso, rainfall data from several meshes will be selected as representative input variables, and rainfall data from meshes that are very similar to the rainfall in the mesh corresponding to the selected input variables will be ignored. Because the rainfall data is very similar, there is usually no problem with selecting rainfall from representative meshes in this way as input variables. However, in cases where there has been frequent localized heavy rainfall in recent years, for example, if there is localized rainfall where the rainfall in the selected mesh happens to be not very heavy, but a nearby mesh that was not selected as an input variable experiences very heavy rainfall, selecting a certain mesh as the representative mesh could result in an erroneous prediction that inflow volume will not increase.

On the other hand, when the input variables described here are integrated, the rainfall in meshes with similar effects is averaged, which means the average rainfall in the integrated meshes. Therefore, even in such cases, it is expected that reasonable forecast results will be obtained, while also improving rational interpretability. The ability to collectively interpret input variables to be integrated in this way, making them easier to explain, is another major advantage of a modular structure (constructive approach).

Finally, the predicted output calculated (synthesized) by the composite output variable prediction unit 7 is sent to the output prediction result observation unit 10 and presented to the user, i.e., the plant manager or operator, as time-series data known as a trend graph. At this time, it is preferable that the current real-time monitored values be displayed simultaneously with the actual output value (stormwater inflow volume) and its predicted value, and that future predicted values be displayed only as the predicted output.

In addition to the synthesized forecast output (forecasted stormwater inflow), the location parameter μy and scale parameter σy can be calculated using the p forecast outputs of the pairwise forecast model, and the range of μy±Kσy (K is a setting parameter set to about 2-3) can be simultaneously displayed, or the maximum value (or the forecast value closest to the 95th percentile, which is close to the maximum value, considering robustness) and minimum value (or the forecast value closest to the 5th percentile, which is close to the minimum value, considering robustness) from the forecast outputs of the pairwise forecast model can be simultaneously displayed as the range of the direct forecast output.

The modular time series data forecasting device of this embodiment introduces a mechanism that enables a rational interpretation of the internal parameters of what are generally called black box models that is consistent with the laws of physics. This improves the model's interpretability and rational explainability, making it closer to white box models, and significantly reduces the complexity of model maintenance, such as partial adjustments, deletions, and integrations.

The fundamental idea behind achieving this effect is to construct a prediction model constructively, in a situation where only input and output data are given and the internal structure is completely unknown, by directly associating elements that may be related to the output with each element (corresponding to a pairwise prediction model) and defining the relationship between each element and the output as the final internal structure of the input and output relationships. By using this idea, the modular time series data prediction device of this embodiment achieves the same explainability and ease of adjustment as a white-box model, without delving into the internal structure, which would otherwise be unknowable.

This method of constructing a modular time series prediction model is well known in the field of statistical machine learning. It is similar to ensemble learning techniques such as bagging and boosting, mentioned above. Bagging and boosting are methods for constructing a highly accurate prediction model called a strong learner (strong learning model) by collecting many weak learners (weak learning models), which are easily constructed but not necessarily highly accurate. Bagging is a method of performing strong learning on weak learners using processes such as majority voting, while boosting is a method of training a strong learner using the results of weak learners. However, while bagging and boosting generate many weak learners by changing the data used for learning (identification), this embodiment differs in that pairwise prediction models equivalent to weak learners are identified by changing the explanatory variables rather than changing the identification data. Note that in the modular time series data prediction device of this embodiment, when defining the prediction method for the composite output, methods such as robust location parameter estimation and weighted averaging using pre-specified weights are similar to the bagging approach, while methods for determining the weights for the weighted average by re-identification are similar to the boosting approach.

Below, several examples of pairwise prediction models in the modular time series prediction device of this embodiment will be described.
(First Example)
FIG. 5 is a diagram illustrating a first example of a pairwise prediction model in a modular time series prediction device according to an embodiment.
FIG. 5 shows an example in which the FIR model or AR model of equation (1) is used as the pairwise prediction model.

The pairwise prediction model MA of this embodiment comprises regression blocks A1-AN and a synthesis block A5 of the predicted output.
For example, each of the regression blocks A1-A3 outputs a weighted moving average component (FIR component/MA component) as a predicted value for the input variable. For example, the regression block AN outputs an autoregressive component (AR component) as a predicted value for the input variable (past output variable).

The prediction output synthesis block A5 acquires the prediction values output from the regression blocks A1-AN, and calculates and outputs a synthesized prediction value (Y(t+1) = F(Y1(t+1), Y2(t+1), ..., Yn(t+1))) by synthesizing the prediction values.
For example, when parameters are identified using a normal ARX model or the like, the weight of the autoregressive component (AR component) becomes large, and if the value of the input variable corresponding to the FIR component changes suddenly, the predicted output may not be able to follow the measured output, resulting in a phase lag in the prediction. However, by adopting the configuration of this embodiment, it is possible to adjust the influence of the AR component, improve the phase lag, and predict changes in the output before the value of the measured output changes.

(Second Example)
FIG. 6 is a diagram illustrating a second example of a pairwise prediction model in a modular time series prediction device according to an embodiment.
FIG. 6 shows a prediction model for intentionally overpredicting and underpredicting, as an example of a prediction model used when it is desired to impart a bias to the combined prediction.

The pairwise prediction model MB of this embodiment includes regression blocks B1-BN and a prediction output synthesis block B5. The synthesis block B5 outputs an overestimated prediction value (Y ₊ (t+1) = F ₊ (Y1(t+1), Y2(t+1), ..., Yn(t+1))), a normal estimation prediction value (Y(t+1) = F(Y1(t+1), Y2(t+1), ..., Yn(t+1))), and an underestimated prediction value (Y _- (t+1) = F _- (Y1(t+1), Y2(t+1), ..., Yn(t+1))).

In this embodiment, the location parameters and scale parameters of the pairwise prediction model are used to make overpredictions and underpredictions in addition to normal predictions, allowing for fine-tuning of prediction results according to needs (requests), such as conservative predictions (predictions that avoid risk).

7 and 8 are diagrams for explaining the effect of the pairwise prediction model of the second embodiment.
FIG. 7 shows an example of the rainwater inflow predicted value and the measured value obtained by Lasso regression with L1 regularization, which is a type of regression method.
FIG. 8 shows an example of the rainwater inflow predicted value and the measured value according to the prediction model of the second embodiment.

Comparing the predicted rainwater inflow values shown in FIG. 7 with the measured values, it can be seen that although the prediction error is small, the estimated values are slightly smaller than the measured values, resulting in an underestimation.
On the other hand, when comparing the results of normal estimation (average estimation), overestimation, and underestimation obtained by the method of this embodiment shown in Figure 8 with the measured values, the normal estimates (predicted inflow volumes) have small prediction errors and good prediction accuracy near peaks, but the predicted values at times when inflow volumes increase are lower than the measured values.

In such cases, if this predicted value is to be used for control, it may be desirable to overestimate the prediction when there is a sudden increase in inflow. In this case, if an overestimated prediction is made, the overall prediction accuracy in terms of MSE (mean squared error) will deteriorate, but it will be possible to make a prediction that suggests the possibility of a sudden increase in inflow when it occurs. Also, if you want to underestimate the normal predicted value, you can adjust the predicted value to be smaller by making an underestimated prediction.

By making it possible to adjust predicted values using such overestimation and underestimation, if, as a result of long-term operation, a tendency (bias) is observed in which normal predicted values tend to underestimate or overestimate actual values, the normal predicted values can be replaced with overestimated or underestimated values to correct the bias and make more accurate predictions.

(Third Example)
FIG. 9 is a diagram illustrating a third example of a pairwise prediction model in a modular time series prediction device according to an embodiment.
In this example, as a pairwise prediction model, an AR model/FIR model such as that in equation (5) is further decomposed, and a prediction model is constructed by decomposing the variables that take time delays into account into a simple regression model in which each variable is regarded as a single explanatory variable.

The pairwise prediction model MC of this embodiment includes prediction output blocks C1-CN, each including multiple time delay blocks, multiple simple regression blocks, and a correlation weighted average block, for each input variable 1-N, and further includes a prediction output value synthesis block C(N+1) that outputs a final synthesis prediction value by synthesizing the prediction values output from the multiple prediction output blocks C1-CN.

The prediction output block C1 includes a plurality (m) of time lag blocks 111-11m, a plurality (m) of simple regression blocks 121-12m, and a correlation weighted average block 13.
The plurality of time delay blocks 111-11m output variables (for example, values at times t-L11, . . . , t-L1m1) delayed by the delay times L11-L1m1 assigned to the input variable 1 (for example, value at time t).

The plurality of simple regression blocks 121-12m use simple regression coefficients to calculate and output predicted values based on the variables output from the time lag blocks 111-11m.
The correlation weighted average block calculates and outputs a correlation weighted average of a plurality of predicted values output from a plurality of simple regression blocks 121-12m.

The predicted output blocks C2-CN have the same configuration as the predicted output block C1, and therefore individual explanations will be omitted.
The synthesis block C(N+1) of the predicted output value synthesizes (for example, performs a correlation weighted average) the predicted values output from the prediction output blocks C1-CN and outputs the final estimated predicted value (Y(t+1) = F(Y1(t+1), Y2(t+1), ..., Yn(t+1))).

The pairwise prediction model of this embodiment can be explained using three terms: delay, gain, and bias, making it easy to fine-tune the parameters. Furthermore, since there is a positive relationship between the simple regression coefficient and the correlation coefficient, the entire prediction model can be constructed simply by calculating the correlation coefficient between each input variable and each output variable, making it easy to interpret the prediction model.
Furthermore, the pairwise prediction model of this embodiment can achieve the same effects as the pairwise prediction model of the first embodiment.

10 and 11 are diagrams for explaining an example of the effect of the pairwise prediction model of the third embodiment.
10 and 11 show an example of rainwater inflow prediction using 10 input variables input to the pairwise prediction model of Example 3. Here, an example of the prediction results for the simplest case is shown, where only one delay time is considered for each input variable.

When estimating the delay time for rainfall intensity (ground rainfall), one of the input variables, an estimated result of 33 minutes was obtained. Figure 11 shows an example of the rainwater inflow forecast results using the above estimation results. The forecast results in Figure 11 show that, overall, predictions were made with a certain degree of accuracy, but there was a large delay in the forecast (phase lag) in areas where the inflow volume increased sharply.

In this case, it can be assumed that the reason for the large phase lag is that the delay time of 33 minutes is too long. Based on this assumption, Figure 10 shows an example of the forecast results for rainwater inflow when the delay time is adjusted to 15 minutes. The forecast results shown in Figure 10 show that by adjusting the delay time, the phase lag during sudden increases in inflow can be significantly improved without significantly affecting the overall forecast accuracy.

In this way, the pairwise prediction model of this embodiment makes it possible to build a prediction model using only easily interpretable parameters, which not only improves interpretability but also makes it easy to fine-tune the prediction model.

It is also known that the simple regression coefficient and correlation coefficient will match if the input and output data are normalized by their respective means and standard deviations; conversely, if they are denormalized, they can be converted into regression coefficients and biases. Therefore, by calculating only the mean, standard deviation, and input/output correlation of each input and output data, it is possible to identify the parameters of all prediction models simply by specifying the delay time.

Furthermore, by calculating the composite output using a weighted average of the correlation coefficient, the only calculations required to build a prediction model are the estimation (specification) of the delay time, the mean and standard deviation of the input and output variables, and the correlation between the input and output variables. All of the parameters included in the prediction model can be explained using only these four estimated values (delay time, mean, standard deviation, and correlation coefficient).

(effect)
As described above, according to this embodiment, when prediction accuracy is insufficient due to reasons such as a change or deterioration in the correlation between input variables and output variables, or the inclusion of a large number of outliers in the time-series data of the input variables, it is possible to easily perform partial adjustments of the pairwise prediction models that affect the synthesized predicted output value or delete parts that have a negative impact, without having to re-identify the entire prediction model. At the same time, this makes it possible to rationally explain the reasons for the deterioration in prediction accuracy.

Furthermore, even when there is extremely strong multicollinearity between input variables and the input variables have roughly the same explanatory power, it is possible to maintain the identifiability of the parameters (unique decidability), making it possible to construct a predictive model that allows for rational explanations that do not contradict the laws of physics.At the same time, by integrating input variables that have roughly the same explanatory power, it is possible to improve interpretability while maintaining prediction accuracy.
That is, according to the present embodiment, it is possible to provide a modular time-series data prediction device, a modular time-series data prediction method, and a program that can rationally explain prediction results and easily adjust prediction results.

Although several embodiments of the present invention have been described, these embodiments are presented as examples and are not intended to limit the scope of the invention. These novel embodiments can be embodied in various other forms, and various omissions, substitutions, and modifications can be made without departing from the spirit of the invention. These embodiments and their modifications are included within the scope and spirit of the invention, and are also included in the scope of the invention and its equivalents as defined in the claims.
(Appendix 1)
An apparatus applied to a system or process having a plurality of process sensors that measure a plurality of process variables at a predetermined cycle,
a data collection and storage unit that includes an output variable data selection unit that collects and stores time series data of the plurality of process variables at a predetermined cycle, and selects at least one output variable to be predicted from the plurality of process variables, and an input variable data selection unit that selects a plurality of input variable candidates from the plurality of process variables;
a pairwise prediction model identification unit that identifies parameters of a pairwise prediction model for each one-input, one-output pair using identification data of the output variable and the plurality of input variables extracted from the time-series data, thereby defining the plurality of pairwise prediction models;
a prediction model synthesis method definition unit that defines a synthesis method for pairwise prediction values output from the plurality of pairwise prediction models;
a pairwise output variable prediction unit that inputs prediction data for the plurality of input variables extracted from the time-series data in a predetermined cycle in the time progression direction or in real time into the plurality of pairwise prediction models, and calculates the pairwise predicted values corresponding to each of the plurality of input variables;
a composite output variable prediction unit that calculates a predicted value of the output variable by combining a plurality of the pairwise predicted values by the combination method;
A modular time series data forecasting device having the above.
(Appendix 2)
2. The modular time series data prediction device according to claim 1, wherein the pairwise prediction model identification unit applies a finite impulse response model to identification data of the plurality of input variables other than the output variable, and applies an autoregressive model to the output variable.
(Appendix 3)
the pairwise prediction model identification unit constructs a one-input, one-output transfer function model using a simple regression model incorporating one or more delay times for each of the plurality of input variables, and defines a model obtained by combining the simple regression models as the pairwise prediction model.
(Appendix 4)
2. The modular time-series data prediction device according to claim 1, wherein the pairwise prediction model identification unit performs nonlinear transformation on identification data for each of the plurality of input variables to construct the pairwise prediction model.
(Appendix 5)
the prediction model synthesis method definition unit defines, as the synthesis method, a weighted average of the plurality of pairwise prediction values, and determines weights in the synthesis method based on an index representing the prediction accuracy of each of the plurality of pairwise prediction models identified by the pairwise prediction model identification unit.
(Appendix 6)
2. The modular time-series data prediction device according to claim 1, wherein the prediction model synthesis method definition unit defines the synthesis method as calculating a weighted average of the pairwise prediction values, and identifies and determines weight values in the synthesis method using the same data as identification data used by the pairwise prediction model identification unit.
(Appendix 7)
a prediction error evaluation unit that evaluates an error between the plurality of pairwise predicted values output from the pairwise output variable prediction unit and the predicted value of the output variable output from the composite output variable prediction unit and an actual value extracted from the time series data at a predetermined cycle or for a predetermined period;
a pairwise prediction model correction unit that corrects the pairwise prediction models used in the pairwise output variable prediction unit by comparing errors in the pairwise prediction values with errors in the predicted values of the output variables to determine the pairwise prediction models to be adjusted, the pairwise prediction models to be integrated, and the pairwise prediction models to be separated;
2. The modular time series data prediction device according to claim 1, further comprising:
(Appendix 8)
and wherein, when at least one of the plurality of errors evaluated by the prediction error evaluation unit exceeds a predetermined allowable error, the pairwise prediction model correction unit re-identifies the pairwise prediction model for the candidates of the input variables corresponding to the error, and corrects the plurality of pairwise prediction models used in the pairwise output variable prediction unit.
(Appendix 9)
8. The modular time-series data prediction device according to claim 7, wherein, when at least one of the plurality of errors evaluated by the prediction error evaluation unit exceeds a predetermined allowable error, the pairwise prediction model correction unit deletes the pairwise prediction model for the input variable corresponding to the error, and corrects the plurality of pairwise prediction models used in the pairwise output variable prediction unit.
(Appendix 10)
the pairwise prediction model correction unit compares the errors evaluated by the prediction error evaluation unit, and, when a difference between the errors is equal to or less than a predetermined value, defines a sum or average of the input variables corresponding to the errors as a new input variable, and corrects the pairwise prediction models so that the pairwise output variable prediction unit uses the pairwise prediction model related to the new input variable.
(Appendix 11)
1. A method for application to a system or process having a plurality of process sensors measuring a plurality of process variables at a predetermined period, comprising:
collecting and storing time series data of the plurality of process variables at a predetermined cycle;
selecting at least one output variable to be predicted from the plurality of process variables;
selecting a plurality of candidate input variables from among the plurality of process variables;
identifying parameters of a pairwise prediction model for each one-input, one-output pair using identification data of the output variable and the plurality of input variables extracted from the time-series data, thereby defining a plurality of the pairwise prediction models;
defining a method for combining pairwise prediction values output from a plurality of said pairwise prediction models;
inputting prediction data of the plurality of input variables extracted from the time-series data in a predetermined cycle or in real time in the time progression direction into the plurality of pairwise prediction models, and calculating the pairwise predicted values corresponding to each of the plurality of input variables;
a combination of the pairwise predicted values by the combination method to calculate a predicted value of the output variable;
Modular time series data forecasting method.
(Appendix 12)
A modular time series data forecasting program that causes a computer to execute the method described in Appendix 11.

1...Urban stormwater drainage process, 11...Inflow meter, 12...Main flow meter, 131-13K...Main water level gauge, 141-14M...Ground rain gauge, 15, 1511-15QP...Radar rain gauge, 17...Inflow culvert, 110...Inflow gate, 18...Rainwater pump well, 19...Rainwater pump, 161, 162...Rainwater pump well water level gauge, 2...Data collection and storage unit, 21...Output variable data selection unit, 22...Input variable data selection unit, 3...Data extraction unit, 31...Data extraction unit for offline prediction model identification, 32...Data extraction unit for online prediction, 33...For evaluation Data extraction unit, 4...pairwise prediction model identification unit, 5...prediction model synthesis method definition unit, 6...pairwise output variable prediction unit, 7...synthetic output variable prediction unit, 8...prediction error evaluation unit, 9...pairwise prediction model correction unit, 10...output prediction result observation unit, 111-11m...time delay block, 121-12m...simple regression block, 13...correlation weighted average block, A1-AN...regression block, A5...synthesis block, B1-BN...regression block, B5...synthesis block, C1-CN...prediction output block, C(N+1)...synthesis block.

Claims (12)

An apparatus applied to a system or process having a plurality of process sensors that measure a plurality of process variables at a predetermined cycle,
a data collection and storage unit that includes an output variable data selection unit that collects and stores time series data of the plurality of process variables at a predetermined cycle, and selects at least one output variable to be predicted from the plurality of process variables, and an input variable data selection unit that selects a plurality of input variables from the plurality of process variables;
a pairwise prediction model identification unit that identifies parameters of a pairwise prediction model for each one-input, one-output pair using identification data of the output variable and the plurality of input variables extracted from the time-series data, thereby defining the plurality of pairwise prediction models;
a prediction model synthesis method definition unit that defines a synthesis method for pairwise prediction values output from the plurality of pairwise prediction models;
a pairwise output variable prediction unit that inputs prediction data for the plurality of input variables extracted from the time-series data in a predetermined cycle in the time progression direction or in real time into the plurality of pairwise prediction models, and calculates the pairwise predicted values corresponding to each of the plurality of input variables;
a composite output variable prediction unit that calculates a predicted value of the output variable by combining a plurality of the pairwise predicted values by the combination method;
A modular time series data forecasting device having the above. The modular time series data prediction device of claim 1, wherein the pairwise prediction model identification unit applies a finite impulse response model to identification data for the plurality of input variables other than the output variable, and applies an autoregressive model to the output variable. 2. The modular time series data prediction device according to claim 1, wherein the pairwise prediction model identification unit constructs a one-input, one-output transfer function model using a simple regression model incorporating one or more delay times for each of the plurality of input variables, and defines a model obtained by combining the simple regression models as the pairwise prediction model. The modular time series data prediction device of claim 1, wherein the pairwise prediction model identification unit performs nonlinear transformation on the identification data for each of the multiple input variables to construct the pairwise prediction model. The modular time series data prediction device of claim 1, wherein the prediction model synthesis method definition unit defines the synthesis method as calculating a weighted average of the plurality of pairwise prediction values, and determines the weights in the synthesis method based on an index representing the prediction accuracy of each of the plurality of pairwise prediction models identified by the pairwise prediction model identification unit. The modular time series data prediction device of claim 1, wherein the prediction model synthesis method definition unit defines the synthesis method as calculating a weighted average of the pairwise prediction values, and identifies and determines the weight values in the synthesis method using the same identification data as used by the pairwise prediction model identification unit. a prediction error evaluation unit that evaluates an error between the plurality of pairwise predicted values output from the pairwise output variable prediction unit and the predicted value of the output variable output from the composite output variable prediction unit and an actual value extracted from the time series data at a predetermined cycle or for a predetermined period;
a pairwise prediction model correction unit that corrects the pairwise prediction models used in the pairwise output variable prediction unit by comparing errors in the pairwise prediction values with errors in the predicted values of the output variables to determine the pairwise prediction models to be adjusted, the pairwise prediction models to be integrated, and the pairwise prediction models to be separated;
2. The modular time series data forecasting device according to claim 1, further comprising: 8. The modular time-series data prediction device according to claim 7, wherein, when at least one of the plurality of errors evaluated by the prediction error evaluation unit exceeds a predetermined allowable error, the pairwise prediction model correction unit re-identifies the pairwise prediction model for the input variable corresponding to the error, and corrects the plurality of pairwise prediction models used in the pairwise output variable prediction unit. A modular time series data prediction device as described in claim 7, wherein, when at least one of the multiple errors evaluated by the prediction error evaluation unit exceeds a predetermined allowable error, the pairwise prediction model correction unit deletes the pairwise prediction model for the input variable corresponding to that error, and corrects the multiple pairwise prediction models used in the pairwise output variable prediction unit. A modular time series data prediction device as described in claim 7, wherein the pairwise prediction model correction unit compares the multiple errors evaluated by the prediction error evaluation unit, and when the difference between the multiple errors is equal to or less than a predetermined value, defines the sum or average of the input variables corresponding to the error as a new input variable, and corrects the multiple pairwise prediction models so that the pairwise output variable prediction unit uses the pairwise prediction model for the new input variable. 1. A method for application to a system or process having a plurality of process sensors measuring a plurality of process variables at a predetermined period, comprising:
The processor:
collecting and storing time series data of the plurality of process variables at a predetermined cycle;
selecting at least one output variable to be predicted from the plurality of process variables;
selecting a plurality of input variables from among said plurality of process variables;
identifying parameters of a pairwise prediction model for each one-input, one-output pair using identification data of the output variable and the plurality of input variables extracted from the time-series data, thereby defining a plurality of the pairwise prediction models;
defining a method for combining pairwise prediction values output from a plurality of said pairwise prediction models;
inputting prediction data of the plurality of input variables extracted from the time-series data in a predetermined cycle or in real time in the time progression direction into the plurality of pairwise prediction models, and calculating the pairwise predicted values corresponding to each of the plurality of input variables;
a combination of the pairwise predicted values by the combination method to calculate a predicted value of the output variable;
Modular time series data forecasting method. A modular time-series data forecasting program that causes a computer to execute the method of claim 11.