WO2011121726A1 - 異常検出装置 - Google Patents

異常検出装置 Download PDF

Info

Publication number
WO2011121726A1
WO2011121726A1 PCT/JP2010/055695 JP2010055695W WO2011121726A1 WO 2011121726 A1 WO2011121726 A1 WO 2011121726A1 JP 2010055695 W JP2010055695 W JP 2010055695W WO 2011121726 A1 WO2011121726 A1 WO 2011121726A1
Authority
WO
WIPO (PCT)
Prior art keywords
value
probability distribution
model
series data
residual
Prior art date
Application number
PCT/JP2010/055695
Other languages
English (en)
French (fr)
Inventor
酢山　明弘
花田　雄一
Original Assignee
株式会社 東芝
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社 東芝 filed Critical 株式会社 東芝
Priority to PCT/JP2010/055695 priority Critical patent/WO2011121726A1/ja
Priority to JP2012507951A priority patent/JP5337909B2/ja
Publication of WO2011121726A1 publication Critical patent/WO2011121726A1/ja
Priority to US13/628,951 priority patent/US8577649B2/en

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0224—Process history based detection method, e.g. whereby history implies the availability of large amounts of data
    • G05B23/0227—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions
    • G05B23/0232—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions based on qualitative trend analysis, e.g. system evolution
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0243—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model
    • G05B23/0254—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model based on a quantitative model, e.g. mathematical relationships between inputs and outputs; functions: observer, Kalman filter, residual calculation, Neural Networks
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B2219/00—Program-control systems
    • G05B2219/20—Pc systems
    • G05B2219/24—Pc safety
    • G05B2219/24042—Signature analysis, compare recorded with current data, if error then alarm

Definitions

  • the present invention relates to an abnormality detection device that detects abnormality of a control device.
  • Model-based anomaly detection is based on the value obtained by modeling the monitored process in advance and taking the residual between the data observed as the process output and the output predicted by the model for a given input. This is a method of outputting the degree of abnormality of the control device.
  • model-based anomaly detection an approximate model with unknown parameters is constructed, and a model learning function for estimating unknown parameters from actual data not including anomalies is included.
  • Patent Document 1 Patent No. 3254624
  • the present invention provides an abnormality detection apparatus capable of detecting an abnormality of a control device having an unknown characteristic with individual differences with high accuracy.
  • the abnormality detection apparatus as one aspect of the present invention is: First to Kth control devices that are arranged at different locations in the medium path for distributing the medium and that respectively output-control the medium in the medium path in accordance with an instruction value given from outside, and the first to Kth An abnormality detection device for first to Kth sensors for measuring an output amount of a medium whose output is controlled by a control device, (Ci + ⁇ i (Xi)) ⁇ (Qi / Qj) ⁇ 2 ⁇ (Cj + ⁇ j (Xj) (Xi is the first observed variable, Bi (Xi) is the first unknown function with Xi as the input variable.
  • Ci is the first constant term
  • Qi is the second observation variable
  • Xj is the third observation variable
  • Bj (Xj) is the second unknown characteristic term that is an unknown function with Xj as the input variable
  • Cj is the first
  • An objective model storage unit for storing an objective model defined by two constant terms, Qj being a fourth observation variable)
  • a first memory that stores first time-series data including the instruction values given to the first to Kth control devices in the first period and measured values by the first to Kth sensors at predetermined time intervals.
  • Xi of the target model is an instruction value of one control device
  • Ci is an unknown constant according to an arrangement location of the one control device.
  • each diagnostic model instance is generated, Define Bi (Xi) and Bj (Xj) in each of the diagnostic model instances generated for each pair as a function that takes a constant value for each section that divides the range of the values of Xi and Xj.
  • the Bi (Xi), the Ci, the Bj (Xj), and the Cj so as to minimize the square of the values of all the diagnostic model instances using the first time series data in the first storage unit.
  • a model optimization unit that obtains each optimized diagnostic model instance by identifying In a second period different from the first period, second time series data including the instruction value given to the first to Kth control devices and a measurement value by the first to Kth sensors is predetermined.
  • a second storage unit for storing at time intervals;
  • a third storage unit for storing accuracy information representing the maximum observation error by the first to Kth sensors and the maximum response error of the first to Kth control devices with respect to the indicated value;
  • For each optimized diagnostic model instance A first residual that is a value of the optimized diagnostic model instance is obtained using the first time-series data in the first storage unit, and a first probability distribution that is a probability distribution of the first residual is generated
  • a second residual that is a value of the optimized diagnostic model instance is obtained using the second time series data in the second storage unit, and a second probability distribution that is a probability distribution of the second residual is generated.
  • the corresponding value in the first time series data is changed according to the accuracy information. Then, using the time series data after the change, obtain a third residual that is the value of the optimized diagnostic model instance, respectively generate a third probability distribution that is a probability distribution of the third residual,
  • the first extended KL (Kulback-Leibler divergence) information amount representing the distance between the second probability distribution and the first probability distribution, and the distance between the second probability distribution and each of the third probability distributions Determining a second extended KL information amount, a determination score calculation unit that calculates a determination score according to a ratio of the first extended KL information amount to a sum of the first and second extended KL information amounts, and For each of the first to Kth control devices and the first to Kth sensors, an overall score that is a minimum value of a determination score obtained from the optimized diagnostic model instance is obtained, and an overall score that is greater than a first threshold is obtained.
  • a comprehensive abnormality determination unit that determines that
  • FIG. 3 is a diagram showing an example different from FIG. It is a figure which shows an example of an input screen. It is a figure which shows an example of an index table. It is a figure which shows an example of a data table. It is a figure which shows the example of several objective models stored in the objective model storage part. It is a flowchart which shows the flow of a process of a model production
  • FIG. 12 is a diagram showing memory contents in the parameter identification calculation process in FIG. It is a figure which shows the mode of a parameter identification calculation. It is a figure which shows the mode of a parameter identification calculation. It is a figure which shows the mode of a parameter identification calculation. It is a figure which shows an example of model information. It is a figure which shows the example of the time series data memorize
  • FIG. 1 is a block diagram schematically showing a configuration of an abnormality detection apparatus as an embodiment of the present invention.
  • FIG. 2 is a diagram for explaining an example for explaining the present embodiment.
  • an air conditioning system that does not have a fan that feeds wind or air (medium) in the duct.
  • VAV VariableVAir Volume
  • CAV AVConstant Air Volume
  • VAV1 to VAV3 are installed.
  • VAV1 to VAV3 are installed in different places.
  • VAV1 to VAV3 are provided with valves X1 to X3 for adjusting the amount of air flowing through each branch duct.
  • the opening instruction value is generated by a controller (not shown) and is given to VAV1 to VAV3.
  • air volume sensors AF1 to AF3 are arranged near VAV1 to VAV3.
  • the air volume sensor observes the air volume (output volume) whose output is controlled by the VAV (control device).
  • the duct configuration is not limited to that shown in FIG. 2, and various modes are conceivable.
  • the number of branch ducts may be an arbitrary number n of two or more.
  • the outlets of the respective branch ducts may exist in a plurality of rooms. In the example of FIG. 3B, each room is adjacent by a releasable door or window.
  • FIG. 1 describes the abnormality detection apparatus in detail.
  • the abnormality detection device 1 includes a model generation unit 11, a model optimization unit 12, an abnormality detection unit 13, a target model storage unit 16, a reference time series data storage unit (first storage unit) 17, a model storage unit 18, and time series data.
  • a storage unit (second storage unit) 19 a result storage unit 20, a sensor information storage unit (third storage unit) 21, a condition input unit 15, and a display unit 16 are provided.
  • Elements 11 to 13 are realized, for example, by causing a CPU to execute a program module.
  • the CPU reads the program module from a program storage unit (not shown), expands it in the memory, and executes it.
  • the program storage unit can be realized by an arbitrary storage medium such as a hard disk device, a memory device, a CD-R, a DVD-R, a DVD-RW, or a flexible disk.
  • Elements 16 to 21 are constituted by an arbitrary storage medium such as a hard disk device, a memory device, a CD-R, or a DVD-ROM.
  • the storage medium may be provided inside the abnormality detection device 1 or may be an external device that can be detached.
  • the condition input unit 15 is an input interface for inputting data by the user, and can be configured by any input means such as a keyboard and a mouse.
  • the display unit 16 is an output interface with the user, and can be configured by an arbitrary display device such as a liquid crystal display device, a plasma display device, an organic EL display device, or a CRT device.
  • the display unit 16 displays data specified by the model generation unit 11 or the abnormality detection unit 13, for example.
  • the condition input unit 15 may include the touch panel function.
  • FIG. 4 shows an example of an input screen for performing input by the condition input unit 15.
  • rooms A to D are selected as building X and floors 2 through 14 as floors. One floor may be selected and one room may be selected.
  • parameters are set for the sensors and actuators to be detected in the room.
  • common settings are made for the rooms A to D on the 2nd to 14th floors.
  • the point object name includes six fields separated by a delimiter (here “.”).
  • the building ID is the ID of the selected building.
  • the floor ID is an ID that identifies the selected floor and room.
  • the large point type is a place type in the room. For example, “I” (Interior) represents the interior of the room, and “P” (Perimeter) represents the window.
  • the point small type is a device type. For example, there are VAV and CAV.
  • the point ID is the ID of the device specified by the small point type.
  • the state quantity indicates a physical quantity or an instruction value observed for the designated device.
  • the reference time series data storage unit (first storage unit) 17 in FIG. 1 stores an index table and a data table in advance.
  • FIG. 5 shows an example of the index table 171.
  • the index table 171 includes a handle ID and a point object name.
  • the index table 171 stores, for example, all point object names and handle IDs that can be set on the input screen of FIG.
  • FIG. 6 shows an example of the data table 172.
  • the data table 172 stores, as reference time series data (first time series data), a history of observation values of the sensor / actuator specified by the corresponding point object name for each handle ID of the index table 171.
  • the reference time series data (first time series data) stored in the data table 171 are all acquired in the first period when the sensor / actuator is normal.
  • the observation values corresponding to each handle ID are simultaneously recorded at regular time intervals (in this example, every 10 minutes).
  • the objective model storage unit 16 stores a plurality of objective models to which the following formula 1 can be applied in the building air conditioning system.
  • the objective model is applicable to a partial structure in the diagnosis target system (in this example, between any two control devices (or between branch ducts) in the building air conditioning system).
  • Equation 1 (Ci + ⁇ i (Xi)) ⁇ (Qi / Qj) ⁇ 2 ⁇ (Cj + ⁇ j (Xj)) included in Equation 1 corresponds to the objective model of the present invention.
  • ⁇ Means power.
  • Formula 1 indicates that Ci, Cj, ⁇ i, ⁇ j of the objective model is obtained so as to minimize the square (R) of the value (e) of the objective model.
  • Equation 1a f (Y) is an operation term for the second and fourth observation variables, and is a function of a known form.
  • f (Y) may be described as f (Qi, Qj).
  • (Ci + ⁇ i (Xi)]) ⁇ f (Y) ⁇ (Cj + ⁇ j (Xj)) corresponds to the target model.
  • (Formula 1a) min_ ⁇ Ci, Cj, ⁇ i, ⁇ j ⁇ R ⁇ (Ci + ⁇ i (Xi)]) ⁇ f (Y)-(Cj + ⁇ j (Xj)) ⁇ ⁇ 2
  • Equation 1 The objective model included in Equation 1 above models the energy conservation law of a medium (fluid or gas) flowing in any two tributary ducts. This energy conservation law is shown in Equation 1b below.
  • Equation 1b The derivation process of the objective model will be briefly described with an example of a duct configuration in which three VAVs are arranged as shown in FIG.
  • Ki1, Ki2, Kj1, and Kj2 are constants representing the pressure loss coefficient of the duct installed in the intelligent / perimeter area of the target zone.
  • I and j indicate two different tributaries and are IDs of ducts (or VAVs).
  • Ki1 and Kj1 are pressure loss coefficients that are steadily lost in the pipelines from the branch points to the respective branch ends to the respective branch ends (the pressure loss coefficients inherent to the branch ducts).
  • Ki2 and Kj2 indicate the pressure loss coefficient of each tributary duct that is variably lost according to the VAV opening of each tributary duct.
  • Qi and Qj are air flow measurement values (m 3 / s) indicated by the air flow sensor.
  • Vi and vj are the valve opening indication values for adjusting the duct air volume (no real unit from 0 to 1).
  • Formula 1b means that the VAV air volume of the tributary j can be estimated from the VAV instruction value (VAV opening) and the VAV air volume of the tributary i and the VAV instruction value (VAV opening) of the tributary j. This relational expression holds true even when the number of VAVs is changed, and is highly versatile because it is less susceptible to disturbances.
  • Equation 1b is derived from Bernoulli's theorem for any two tributary ducts.
  • the derivation process of Equation 1b based on Bernoulli's theorem is omitted for simplicity of explanation.
  • the above generalized objective model can be derived by transforming Equation 1b.
  • FIG. 7 shows a plurality of example objective models stored in the objective model storage unit 16.
  • the character string described in ⁇ > included in each objective model that stores the objective models (entries) of F-01, F02, F03, F04, F05, F06 .... is an abstract representation for the point object name. Means.
  • [S1], [s2], [s3], [s4] and [s5] variables are associated with arbitrary character strings. However, if there are multiple variables in the same entry, each variable must be associated with the same string. For example, the objective model of F-01 includes three [s1], but each of these [s1] must be associated with the same character string.
  • the model generation unit 11 includes an input information analysis unit 111 and an instance generation unit 112. Using these elements 111 and 112, the model generation unit 11 generates an instance of Formula 1 based on the information (point object name) input by the condition input unit 15 and the target model storage unit 16.
  • Xi, Xj, Qi, Qj in Equation 1.
  • Xi, Xj, Qi, Qj in Equation 1.
  • Opening Indication Value is assigned to Xi in Formula 1
  • X Building.2A.I.VAV.1 Air Volume is assigned to Qi.
  • FIG. 8 is a flowchart showing a process flow of the model generation unit 11.
  • the model generation unit 11 starts processing with a user input of information from the condition input unit 15 as a trigger.
  • the input information reading unit 111 reads the information (a plurality of point object names) registered by the condition input unit 15 and the index table 171 of the reference time series data storage unit 17 into the memory (S11).
  • the top three information (building ID, floor ID, point major classification ID) in each point object name is extracted, and a list of elements (partial configuration list) including the extracted top three information is created (S12).
  • step S13 an element that has not yet been selected is selected from the partial configuration list.
  • step S14 a diagnostic model instance is generated for the element selected in step S13 (S14).
  • the selected element is “X Building 2A.I”.
  • [s1] is “X building”
  • [s2] is “2A”
  • [s3] is uniquely determined by “I”
  • [S4] and [s5] both correspond to point IDs.
  • ([s4], [s5]) (1, 2), (2, 1), (1, 3), (3, 1), (2, 3), (3, 2) Therefore, the following six instances (diagnostic model instances) ⁇ 1> to ⁇ 6> are obtained for the model of F-01 (see FIG. 9).
  • step S15 corresponding to each diagnostic model instance generated in step S14, a partial diagnosis list (parameter table and access table) for matching with the target model is registered.
  • the parameter table is a table storing the coefficients (Ci, Cj, ⁇ i, ⁇ j) of each diagnosis model instance.
  • An example of a parameter table is shown on the right side of FIG. 9 (linked to the “PARMS” attribute).
  • the coefficient C [3] is the coefficient C1 for F-01
  • the coefficient C [2] is the coefficient C2 for F-01
  • the coefficient is the coefficient ⁇ 1 for F-01 ⁇ [2] is stored as the coefficient ⁇ 2 of ⁇ [3] and F-01. That is, for the bottom instance, “C [3] 2, C [2” corresponding to C1, C2, ⁇ 1, ⁇ 2 of F-01 (or Ci, Cj, ⁇ 2, ⁇ 2 in Formula 1) ], ⁇ [3], ⁇ [2] ”are registered in the parameter table. At this time, the values of these parameters have not been identified (unknown).
  • the access table stores IDs corresponding to observation variables (opening instruction value, air volume) included in each diagnosis model instance.
  • An example of an access table is shown on the right side of FIG. 9 (linked to the “ACCESS” attribute).
  • the access table is used for mapping to reference time-series data stored in the reference time-series data storage unit 17 and mapping to time-series data stored in the time-series data storage unit 19 used in abnormality diagnosis described later. It is done. For example, for instance ⁇ 6> ⁇ X Building 12C.I.VAV.3 Opening Indication Value> ⁇ X building 12C.I.VAV.2.Opening indication value> ⁇ X Building 12C.I.VAV.3.Airflow> ⁇ X Building 12C.I.VAV.3.Airflow> There are four observation variables The corresponding handle IDs are 135, 134, 145, and 144 (see the index table 171 in FIG. 5). Therefore, “135,134,145,144” is registered in the access table for the instance ⁇ 6>.
  • steps S14 and S15 described above are performed in the same manner for all the remaining elements of the partial configuration list.
  • a model list including model information for each floor of the building and each floor is stored in the model storage unit 18 (step S17).
  • the model list is the output of the model generator and is shown in FIG.
  • FIG. 9 shows an example of model information generated for “12C” (room C on the 12th floor). Similar model information exists for the remaining floors and rooms. A set of these model information corresponds to a model list.
  • the model list includes model information of each floor and each room selected in the condition input unit, and the numbers from 1 to Nm are model information numbers. A number is set for each set of building ID and floor ID, and the number from 1 to Nm increases sequentially as model information is generated.
  • FIG. 10 shows an example of model information generated for “2A” (room A on the second floor). Since the model information “2A” is created first in the above processing flow, number 1 is assigned, and no other model information exists yet at this point. Generation of model information for “12C” is performed 39th, and number 39 is assigned (see FIG. 9). In this example, the last number Nm is 52.
  • each model information includes the following attributes.
  • the model information number of the room C on the 12th floor of the X building is 39, and “ID” indicates a link to the identifier “X building.12C” of the model information.
  • Fig indicates a link to the file in which the screen data of the input screen by the condition input unit 15 is stored.
  • “Bind” indicates a link to the diagnosis model instance described above.
  • Attribute “C” Constant term
  • “segment” segment
  • “B” unknown characteristic term
  • “w” confidence value
  • the model optimization unit 12 in FIG. 1 includes a data reading unit 121, an unknown characteristic term initial calculation unit 122, a constant term convergence calculation unit 123, and an unknown characteristic term convergence calculation unit 124.
  • the model optimization unit 12 minimizes the square of the value of the diagnostic model instance (R in Formula 1) with respect to the reference time series data stored in the reference time series data storage unit 17.
  • the parameter of the diagnostic model instance is learned.
  • the data reading unit 121 reads the reference time series data stored in the reference time series data storage unit 17 into the memory.
  • the unknown characteristic term initial calculation unit 122 fixes the constant term included in each diagnostic model instance and makes an error (the diagnostic model instance).
  • the unknown characteristic term is identified and stored in the memory so that the square of (value) is minimized.
  • the constant term convergence calculation unit 123 fixes the unknown characteristic term in the state currently stored in the memory, and identifies the constant term so that the square of the error is minimized with respect to the reference time series data read into the memory. Store in memory.
  • the unknown characteristic term convergence calculation unit 124 fixes the constant term in the state currently stored in the memory and sets the unknown characteristic term so that the square of the error is minimized with respect to the reference time series data read into the memory. Identify and store in memory.
  • the constant term convergence calculation unit 123 and the unknown characteristic term convergence calculation unit 124 repeatedly execute processing alternately until the end condition is satisfied, and store the parameters at that time in the model storage unit 18 when the end condition is satisfied.
  • the diagnostic model instance in which the parameters thus obtained are set is referred to as an optimized diagnostic model instance.
  • FIG. 11 is a flowchart showing the flow of processing by the model optimization unit 12.
  • the model optimization unit 12 is executed after the model generation unit 11. Therefore, before the model optimization unit 12 is executed, model information for the floor / room selected by the user of the X building is stored in the model storage unit 18.
  • model information of an arbitrary number in the model list stored in the model storage unit 18 (S21). For example, the model information of No. 39 “X building, 12th floor, room C” is read.
  • the corresponding data from the reference time series data in the reference time series data storage unit 17 (the observed value corresponding to the handle ID of the ACCESS attribute and the time The stamp (time information) is read into the memory (S22).
  • Missing value processing is performed by an existing method. For example, it is desirable to apply a modified median filter that does not make a large value change to input data, has a relatively high ability to remove missing values and outliers, and has a small amount of calculation.
  • steps S24 to S28 are performed by the unknown characteristic term initial calculation unit 122, the constant term convergence calculation unit 123, and the unknown characteristic term convergence calculation unit 124.
  • Step S24 is performed by the unknown characteristic term initial calculation unit 122.
  • ⁇ k (X) which is an unknown characteristic term, is not given by a predetermined function
  • ⁇ k (X) is replaced with ⁇ [k, s (X)].
  • ⁇ [k, s (X)] is a function that gives a constant for each of the M segments obtained by dividing the range of possible values of X by a certain width.
  • s (X) is a function that returns the index of the segment to which X belongs when X is given, and outputs any natural number from 1 to M for the given X (M: segment number).
  • k is an unknown characteristic term identifier corresponding to i or j in the example of Equation 1.
  • ⁇ [k, s (X)] when given X, returns the constant of the segment to which X belongs.
  • Formula 1 that is an expression indicating optimization can be expressed by the following Formula 2.
  • Equation 2 is C [i], C [j], ⁇ [i, 1],..., ⁇ [i, M], ⁇ [j, 1],..., ⁇ [j, M] ⁇ ⁇ ⁇ means to obtain (optimize) a total of (2 ⁇ M + 2) parameters.
  • C [i]> 0 and C [j]> 0 are constraints.
  • C1, C2, C3 can be said to be unknown constant terms corresponding to the locations (ducts) where VAV1 to 3 are arranged.
  • ⁇ [1,1: M], ⁇ [2,1: M], ⁇ [3,1: M] corresponding to the duct pressure loss coefficient estimated from the opening indication values of VAV1 to VAV3
  • ⁇ [0,1: M] that is, ⁇ [0,1],..., ⁇ [0, M]. That is, ⁇ [1,1], ⁇ [2,1], ⁇ [3,1] all have the same value ⁇ [0,1], and ⁇ [1,2], ⁇ [2,2], ⁇ [3,2] all have the same value ⁇ [0,2],... ⁇ [1, M], ⁇ [2, M], ⁇ [3, M] all have the same value ⁇ [ Suppose it has 0, M].
  • ⁇ [0,1: M] corresponding to the common tributary duct pressure loss coefficient is obtained so that the square error with respect to the data read from the reference time series data storage unit 17 is minimized.
  • FIG. 13 shows the calculation in step 1 (S25).
  • the X axis is the VAV opening instruction value
  • the Y axis is ⁇ [0, s (xk)] + C0.
  • Fig. 12 shows the memory contents in the calculation process of steps 1 to 3 (S25 to S27), and variables to be optimized are held as an array table on the memory.
  • the memory contents of this step 1 are shown in the upper part of FIG. While the variables of the cells that are outlined are fixed, the variables of the cells that are shaded are optimized.
  • Equation 4.1 is an instance of Equation 4, (S26).
  • C1, C2, and C3 can be found by solving (S26).
  • FIG. 14 shows the state of the calculation in step 2 (S26).
  • C1 old , C2 old , C3 old on the left of the figure are C1, C2, C3 before optimization by this step
  • C1 new , C2 new , C3 new on the right of the figure are C1, old after optimization by this step C2, C3.
  • the X axis is a VAV opening instruction value, specifically, x1 is a VAV1 opening instruction value, x2 is a VAV2 opening instruction value, and x3 is a VAV3 opening instruction value.
  • the Y axis is ⁇ [k, s (x)] + Ck.
  • Step 2 The memory contents in this step 2 (S26) are shown second from the top in FIG. While the variables of the cells that are outlined are fixed, the variables of the cells that are shaded are optimized. In Step 2, C1, C2, and C3 are optimized while ⁇ [1,1: M], ⁇ [2,1: M], and ⁇ [3,1: M] are fixed as described above.
  • Step 3 the unknown characteristic term convergence calculation unit 124 finds the following Expression 5.1, which is an instance of Expression 5.
  • Fig. 15 shows the calculation in step 3 (S27).
  • ⁇ old [1, s (x1)], ⁇ old [2, s (x2)], ⁇ old [3, s (x3)] on the left of the figure are ⁇ [1, s ( x1)], ⁇ [2, s (x2)], ⁇ [3, s (x3)], ⁇ new [1, s (x1)], ⁇ new [2, s (x2)], ⁇ new [3, s (x3)] are ⁇ [1, s (x1)], ⁇ [2, s (x2)], and ⁇ [3, s (x3)] after optimization in this step.
  • Step 3 The memory contents in this step 3 (S27) are shown second from the bottom in FIG. While the variables of the cells that are outlined are fixed, the variables of the cells that are shaded are optimized. In Step 3, ⁇ [1,1: M], ⁇ [2,1: M], and ⁇ [3,1: M] are optimized while C1, C2, and C3 are fixed as described above.
  • step S27 the following processing may be performed to improve the reliability of the abnormality determination described later.
  • the partial decomposition ⁇ [k, s] obtained from the sparse part data is a value to be obtained truly There is a possibility that the error becomes large.
  • a confidence value w [k, s] for ⁇ [k, s] is set so that processing based on the reliability of ⁇ can be performed in the subsequent abnormality determination.
  • the greater the number of data included in the section s ( 1, 2,..., M), the higher the reliability value.
  • a function or table that provides a higher confidence value as the number of data increases is prepared, and the confidence value of the section is calculated by giving the function or table the number of data included in the section.
  • step S27 The memory contents when the confidence value is calculated in step S27 are shown in the lower part of FIG. Since the calculation of the confidence value only needs to be performed once, it is not necessary to repeat it every time step S27 is executed.
  • Ck 1, 2, 3
  • the termination condition is determined according to the solver used. For example, “convergence to solution”, “variation amount of objective variable is smaller than allowable range”, “number of iterations reaches maximum” Can be mentioned.
  • Fig. 16 shows an example of model information after the parameter information is added in step S29.
  • the time-series data storage unit (second time-series data storage unit) 19 sequentially stores the state quantities (instruction value, air volume, etc.) of the devices subject to abnormality detection at time intervals (second time-series data). Data) is memorized in the lil time. Specifically, the state quantity corresponding to each handle ID as shown in FIG. 5 is sequentially stored.
  • An example of time-series data stored in the time-series data storage unit 19 is shown in FIG.
  • the period in which the time series data storage unit 19 stores data is a second period different from the first period stored in the reference time series data storage unit 17. The second period may be after or before the first period.
  • the sensor information storage unit 21 stores sensor / actuator accuracy information.
  • the accuracy information for example, the accuracy guaranteed by the manufacturer at factory shipment is used.
  • the maximum detection error for each air volume sensor is stored in the sensor information storage unit 21 as ⁇ 5% or ⁇ 0.5 (m / s), for example.
  • the maximum response error of the actuator with respect to the opening degree instruction value is stored in the sensor information storage unit 21.
  • the maximum response error of the actuator indicates how much the error occurs in the actual operation opening of the valve with respect to the indicated value.
  • the accuracy information of the actuator may be stored in the sensor information storage unit 21 after being converted into an expression with an opening degree instruction value. For example, an opening instruction value range corresponding to the opening degree error range of the valve degree may be obtained, and the obtained range may be stored as accuracy information.
  • the abnormality detection unit 13 includes a model reading unit 131, a data selection unit 132, an abnormality score calculation unit (determination score calculation unit) 133, a general abnormality determination unit 134, and an abnormality detection presentation information generation unit 135.
  • the abnormality detection unit 13 reads time series data for an arbitrary period to be diagnosed from the time series data storage unit 19, and determines whether there is an abnormality in the control device (actuator) / sensor from the read time series data To do. Then, the determination result is stored in the result storage unit 20 and displayed on the display unit 16.
  • model generation unit 11 and the model optimization unit 12 described so far are blocks that are executed immediately after the introduction of the abnormality detection device 1 or when the environment surrounding the device that is the target of the abnormality detection has changed.
  • the abnormality detection unit 13 is executed constantly or at regular intervals (for example, once daily). Before executing the process of the abnormality detection unit 13, it is assumed that the model information generated and updated by the model generation unit 11 and the model optimization unit 12 is stored in the model storage unit 18.
  • FIG. 18 is a flowchart showing the flow of processing by the elements 131 to 134 in the abnormality detection unit 13.
  • the model reading unit 131 reads model information stored in the model storage unit 18 (S301). More specifically, model information (each room on each floor) stored in the model storage unit 18 is read using the building ID as a key. The following steps are performed for each model information, and abnormality diagnosis is performed. For ease of explanation, explanation will be made assuming model information of room C on the 12th floor.
  • the data selection unit 132 reads the corresponding time series data (value and time corresponding to the ACCESS attribute) from the time series data storage unit 19 to the memory using the ACCESS attribute for each diagnosis model instance of the model information ( S303).
  • the data selection unit 132 performs missing value processing on the read time series data (S304).
  • An existing method may be used for the missing value processing. It is desirable to apply a modified median filter that does not make a large change in the read time-series data, has a relatively high ability to remove missing values and outliers, and has a small amount of calculation. This step can be omitted.
  • the abnormality score (determination score) is calculated by the abnormality score calculation unit 133 (S305). This process is repeatedly executed for each diagnostic model instance in the model information.
  • FIG. 19 is a flowchart showing the flow of abnormality score calculation processing by the abnormality score calculation unit 133.
  • Equation 6 residual calculation equation
  • Equation 7 By calculating (reliability calculation formula), a residual vector and a reliability vector are obtained (S305a, S305b). Details will be described below with reference to FIGS. 20 and 21.
  • ID133 is the opening instruction value of VAV1
  • ID 134 is the opening instruction value of VAV2
  • ID 143 is the air volume of VAV1
  • ID 144 is the air volume of VAV2 (see FIG. 5).
  • each unknown characteristic term B1 and B2 is 0.98 corresponding to the segment attribute segment “0.8”.
  • the air volume values Q1 and Q2 of the air volume sensors AF1 and AF2 are 0.20 and 0.21 (time series data values of ID143 and 144, respectively).
  • the reliability of the diagnostic model instance ⁇ 1> is calculated as follows.
  • Equation 6 and Equation 7 shown above are performed in the same manner for the data after the second time in the time series data.
  • the residual vector e ⁇ 1> in which the residuals e for the number of data (time number) included in the time series data are arranged for the diagnostic model instance ⁇ 1> and the same number of reliability r are obtained.
  • the numerical value in ⁇ > indicates that it corresponds to the first diagnostic model instance ⁇ 1>.
  • step S305c the value of each element in the reliability vector is inspected, and all elements having a value equal to or greater than a predetermined reliability threshold are specified. Then, a residual corresponding to the same position (time) as the identified element is extracted from the residual vector. For example, if rx (x is an arbitrary value from 1 to n) is equal to or greater than the confidence threshold, the residual ex is extracted from the residual vector. Residuals that have not been extracted are not used in the subsequent processing because they have low reliability.
  • the confidence threshold corresponds to the second threshold of the present invention.
  • the corresponding accuracy information is read from the sensor information storage unit 21, and an abnormal score (determination score) is obtained using the read accuracy information and the residual group extracted in step S305c (S305d).
  • the accuracy information to be read out is, for example, accuracy information of the air volume sensors AF1 and AF2 and accuracy information of the VAV1 and VAV2 if the diagnosis model instance ⁇ 1>.
  • FIG. 22 is a diagram for explaining a method for calculating an abnormal score.
  • the probability distribution P l (El) of the residual group extracted in step S305c is obtained. More specifically, the mean and variance of the residual e are calculated, and a probability distribution P l (el) is obtained based on the mean and variance assuming a normal distribution. That is, P l is a probability distribution of a residual obtained by using time series data read from the time series data storage unit 19 and using a parameter of the l (el) th diagnostic index.
  • the probability distribution P l (L) corresponds to the second probability distribution of the present invention.
  • a probability distribution P n, l (el) when there is no abnormality is obtained.
  • a residual group is extracted in steps S305b and S305c.
  • the average and variance of the extracted residual group are obtained, and a probability distribution P n, l (el) is obtained based on the average and variance assuming a normal distribution. That is, P n, l is a residual probability distribution obtained by using the reference time series data read from the reference time series data 17 and using the parameters of the l (el) th diagnostic index.
  • the probability distribution P n, l (L) corresponds to the first probability distribution of the present invention.
  • a probability distribution Pa , l (el), k is obtained when it is assumed that an abnormality has occurred statically for each sensor / VAV (actuator).
  • k is the column number of the access table.
  • the probability distributions P a, l (el), 1 , P a, l (el), 2 , P a, l (el), 3 , Pa , l ( El), ask for 4 .
  • the value of the l (el) row k column of the access table is corrected according to the corresponding accuracy information.
  • the probability distribution Pa, l (el), k corresponds to the third probability distribution of the present invention.
  • Fig. 23 shows how the values in the k column of the access table are corrected in the reference time series data.
  • the maximum error in the plus direction (or minus direction) of the accuracy information is added here.
  • the maximum error +0.05 is added.
  • the instruction value is corrected by obtaining a difference between the instruction values corresponding to the maximum error of the opening, and correcting the obtained difference to the instruction value.
  • the correction is performed in the same way in (2) to (4). The case where both the plus direction and the minus direction of accuracy information are considered will be described later.
  • the abnormal score is calculated using the extended KL information (Kulback-Leibler divergence).
  • KL information is used as a measure for the difference between probability distributions.
  • the KL information amount is not a distance because it has no symmetry, but the extended KL information amount used here has a symmetry and can be defined as a distance between probability distributions.
  • the extended KL information amount regarding the probability distributions P and Q of the discrete random variables is defined by Equation 9 below. It is known that the extended KL information amount is 0 when the probability distributions P and Q match, and the value increases as the deviation increases and does not take a negative value.
  • the abnormal score Score (l, k) for the k th sensor / actuator is obtained as follows in the l (optimum) optimized diagnosis model instance.
  • a deviation (expanded KL information amount) between the probability distribution P l (el) and the probability distribution P n, l is obtained by the following formula 10.
  • the probability distribution P l (el) and the deviation of the probability distributions Pa , l (el), k (expanded KL information amount) are expressed by the following equation 11: Ask for.
  • the anomaly score Score (l, k) for the kth sensor can be defined as follows.
  • the probability distribution P n, l (el) corresponding to the reference time series data is also calculated as the square of the residual e (e ⁇ 2), the mean and variance of the e ⁇ 2 group are calculated, and the chi-square distribution is assumed. Then, the probability distribution is obtained.
  • the probability distribution P l (el) corresponding to the time series data in the time series data storage unit 19 also calculates the square (e ⁇ 2) of the residual e, calculates the mean and variance of the e ⁇ 2 group, A probability distribution is obtained assuming a square distribution. Thereafter, in the same manner as before, the amount of extended KL information is obtained, and the abnormal score is calculated.
  • FIG. 24 is a flowchart showing the flow of the comprehensive abnormality determination process by the comprehensive abnormality determination unit 134.
  • a total score is calculated for each sensor / actuator (VAV1 to VAV3, air volume sensors AF1 to AF3) (S306a).
  • the total score is calculated using Equation 13 below.
  • the minimum value among the abnormality scores obtained for each l is used as the total score.
  • Fig. 25 shows an example of calculating the total score.
  • the overall score of the VAV1 opening instruction value (ID133) is the minimum value 0.08 among the abnormal scores 0.96, 0.95, 0.08, and 0.09 calculated by the diagnostic indices ⁇ 1> to ⁇ 4>.
  • the total score is obtained by selecting the minimum value of the abnormal scores.
  • step S306b the total score of each sensor / actuator is compared with a predetermined determination threshold, and if the total score is greater than the determination threshold, an abnormal label is set for the corresponding sensor / actuator (S306c). When it is below the threshold, a normal label is set (S306d).
  • an abnormal label is set for VAV2 with an overall score of 0.92, and normal labels are set for the other VAVs 1, 3 and the air volume sensors AF1 to AF3.
  • the determination threshold corresponds to the first threshold of the present invention.
  • step S306e the diagnosis target period (period using time-series data), the calculated abnormality score, the total score, and the label set for the sensor / actuator are written in the result storage unit 20 as the diagnosis result of the abnormality detection unit 13. (S306e).
  • the abnormality score of the sensor / actuator is obtained for each diagnosis model instance ⁇ 1> to ⁇ 6>, and the overall score is obtained by selecting the smallest abnormality score for each sensor / actuator. Then, an abnormal label is set for the sensor / actuator whose total score is larger than the determination threshold.
  • a normal score (determination score) of the sensor / actuator is obtained for each diagnosis model instance ⁇ 1> to ⁇ 6> by the following formula 12a. The maximum normal score is selected for each sensor / actuator, an abnormal label is set for the sensor / actuator (point object name) whose overall score is smaller than the judgment threshold, and the sensor / actuator whose overall score is larger than the judgment threshold is normal.
  • a label may be set.
  • the abnormality detection presentation information generation unit 135 of the abnormality output unit 13 uses information stored in the time series data storage unit 19, the reference time series data storage unit 17, the result storage unit 20, and the accuracy information storage unit 21, A result display screen is dynamically generated according to a predetermined format, and the generated screen is displayed on the display unit 16.
  • FIG. 26, FIG. 27, and FIG. 28 are screen examples of the result display generated by the abnormality detection presentation information generation unit 136.
  • Fig. 26 shows a display example when "List" in the screen is selected.
  • the point object name and the total score for each VAV are displayed in descending order of the total score. However, in the illustrated example, the display of the building name of the point object name is omitted.
  • the rooms equipped with sensors and VAVs (actuators) with abnormal labels are shaded.
  • Fig. 27 shows a display example when "Layout" on the screen is selected.
  • the screen shown in the figure shows a display example when room C on the 12th floor is selected from the shaded rooms in FIG.
  • the point object name “X building 12.C.I.VAV.2. Opening indication value” related to VAV2 to which an abnormal label is set is shaded.
  • FIG. 28 shows a display example when “Graph” in the screen is selected.
  • Fig. 28 shows a display example when room C on the 12th floor is selected.
  • a coordinate system is shown in which the horizontal axis is time and the vertical axis is an opening instruction value.
  • an opening instruction value graph G101 for VAV1 an opening instruction value graph G201 for VAV2, and an opening instruction value graph G301 for VAV3 are displayed.
  • the horizontal axis shows the time and the vertical axis shows the airflow coordinate system.
  • this coordinate system there are a graph G102 of the airflow value of VAV1 (value of the airflow sensor AF1), a graph G202 of the airflow value of VAV2 (value of the airflow sensor AF2), a graph G302 of the airflow value of VAV3 (value of the airflow sensor AF3) It is displayed.
  • Each graph G101, G201, G301, G102, G202, G302 uses the date information used in the diagnosis (in this example, September 1, 2009, 03:00 to 23:59), and the time series data storage unit 19 Can be obtained by reading the time-series data for this period (in this example, one day) and plotting the read data.
  • Abnormal section Z1 is a section where an abnormal opening of VAV2 is detected.
  • the analysis result screen of FIG. 29 is generated and displayed.
  • This screen is also dynamically generated using the information stored in the time-series data storage unit 19 and the result storage unit 20 by the abnormality detection presentation information generation unit 136, as in FIGS.
  • Opening indication value has an abnormal label (overall score 0.92).
  • model information of X building .12C is read, and a diagnostic index ⁇ 2> (see FIG. 25) that outputs an overall score of 0.92 is referred to.
  • the time-series data (for the period to be diagnosed) corresponding to the ACCESS index numbers 133, 134, 143, and 144 of the diagnosis index ⁇ 2> are read from the time-series data storage unit 19 and the reference time-series data storage unit 17, respectively.
  • Graphs G101 and G201 at the upper left of the screen and G102 and 202 at the lower left are plotted with data (observed values) read from the time series data storage unit 19, and are the same as those shown in FIG.
  • the graph G211 at the top left of the screen is the estimated opening instruction value for VAV2.
  • This estimated opening instruction value is obtained using Equation 6. Specifically, e in Equation 6 is set to 0, and the parameter is substituted using the PARAMS information of the diagnostic model instance. Then, the time series data read from the time series data storage unit 19 is substituted into variables x [j], Q [i], Q [j] other than x [i] determined to be abnormal, and x [i ] x [i] is an estimated value (period data estimated value).
  • the regression plot is shown on the right side of the screen. This regression plot is created as follows. Similarly, x [i] is compared by substituting the reference time series data read from the reference time series data storage unit 17 into the variables x [j], Q [i], and Q [j] in Equation 6. Obtain as an estimate of the data. Then, “ ⁇ ” is plotted at the intersection of the value related to the I.VAV2. Opening degree indication value of the read reference time series data and the estimated value of the comparison data. On the other hand, “x” is plotted at the intersection of the I.VAV2. Opening instruction value of the time series data read from the time series data storage unit 19 and the period data estimated value obtained above.
  • diagnostic model instances ⁇ 1> and ⁇ 2> are generated for a pair of VAV1 and VAV2.
  • FIG. 2 anomaly detection is possible even if only one of them is generated for each pair. Therefore, in this case, only one of ⁇ 1> and ⁇ 2> is generated for the pair of VAV1 and VAV2, and ⁇ 3> and ⁇ 4> for the pair of VAV1 and VAV3. Any one of these may be generated, and only one of ⁇ 5> and ⁇ 6> may be generated for the pair of VAV2 and VAV3.
  • all pairs of two VAVs are created in accordance with the target model F01 and a diagnostic model instance is generated for each pair.
  • the present invention is not limited to this.
  • a diagnostic model instance does not have to be generated for a pair that is known in advance to have a weak mutual influence.
  • a pair for which a diagnostic model instance should be generated may be specified in advance.
  • the objective model F01 implies that all pairs are generated.
  • the opening instruction value of the VAV installed in the duct and the VAV air volume are taken as examples.
  • the present invention is not limited to this, and can be applied to the following examples, for example. It is.
  • Equation (1) an example of the cold water valve opening indication value (vw) and the water amount (Qw) of the cold water valve in the water amount adjusting device is also possible.
  • the following model instance can be generated from Equation (1).

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Abstract

［課題］個体差のある未知の特性を持つ制御機器等の異常検出を高い精度で行う。 ［解決手段］本発明は、第1の期間において第1～第K制御機器に与えられた指示値と、第1～第Kセンサによる計測値とを含む第1時系列データを記憶する第1記憶部と、前記第1～第K制御機器における2つの制御機器のペア毎に所定の目的モデルの診断モデルインスタンスを生成し、前記診断モデルインスタンスのパラメータを同定した最適化された診断モデルインスタンスを得るモデル最適化部と、第2の期間に取得された第2時系列データを記憶する第2記憶部と、前記最適化された診断モデルインスタンス毎に前記第1および第2時系列データを用いてそれぞれ関連する制御機器およびセンサに対する判定スコアを計算する判定スコア計算部と、前記第1～第K制御機器および前記第1～Kセンサ毎に前記判定スコアに基づき異常有無の判定を行う総合異常判定部とを備える。

Description

異常検出装置
　本発明は、制御機器の異常検出を行う異常検出装置に関する。
　制御機器の異常検出の方法の1つとして、モデルベース異常検出がある。モデルベース異常検出とは、予め監視対象プロセスをモデル化しておき、ある入力に対して、プロセスの出力として観測されたデータと、モデルによって予測された出力との残差をとり、その値に基づいて制御機器の異常度合を出力する方法である。モデルベース異常検出では、未知パラメータを持つ近似モデルを構築し、異常を含まない実績データから未知パラメータを推定するモデル学習機能を含んでいる。
特許第3254624号
　モデルベースで異常検出を行う場合、バルブの開度に対するCv特性(流量特性)など、監視対象とするプロセスに用いられる機器特性を、パラメトリックな関数式として定義し、近似モデルを構築する必要が多々ある。
　しかしながら、機器特性は、数式として与えられることは稀であるため、新しい機器が導入される毎に、精度が高い近似モデルを与えることはコスト高である。また、各機器には固体差が存在するため、特性グラフからバイアスがかかった形ではずれる問題も存在する。
　たとえば特許文献1(特許第3254624号)では、調節弁のスティックスリップの異常検出を行うための近似モデルとして、摺動特性を x’’ = A * x’ + B として定義し、パラメータA,Bを実績データから推定している。
　本発明は、個体差のある未知の特性を持つ制御機器の異常検出を高い精度で行うことを可能とした異常検出装置を提供する。
　本発明の一態様としての異常検出装置は、
　媒体を流通させる媒体経路における異なる箇所に配置され、それぞれ外部から与えられた指示値に従って前記媒体経路内の前記媒体を共通に出力制御する第1～第K制御機器と、前記第1～第K制御機器によって出力制御される媒体の出力量を計測する第1～第Kセンサに対する異常検出装置であって、
　(Ci + βi(Xi) ) × (Qi/Qj)^2 － (Cj + βj(Xj) (Xiは第1観測変数、Bi(Xi)は Xiを入力変数とする未知の関数である第1未知特性項、Ciは第1定数項、Qiは第2観測変数、Xjは第3観測変数、Bj(Xj)は Xjを入力変数とする未知の関数である第2未知特性項、Cjは第2定数項、Qjは第4観測変数)によって定義される目的モデルを記憶する目的モデル格納部と、
　第1の期間において前記第1～第K制御機器に与えられた前記指示値と、前記第1～第Kセンサによる計測値とを含む第1時系列データを所定時間間隔で記憶する第1記憶部と、
　前記第1～第K制御機器における2つの制御機器のペア毎に、前記目的モデルの前記Xiに一方の制御機器の指示値、前記Ciに前記一方の制御機器の配置箇所に応じた未知の定数、前記Qiに前記一方の制御機器に対応するセンサの計測値、前記Xjに他方の制御機器の指示値、前記Cjに前記他方の制御機器の配置箇所に応じた未知の定数、前記Qjに前記他方の制御機器に対応するセンサの計測値を割り当てることによって、それぞれ診断モデルインスタンスを生成し、
　前記ペア毎に生成したすべての診断モデルインスタンスのそれぞれにおける前記Bi(Xi)およびBj(Xj)を、前記XiおよびXjの値の範囲をそれぞれ分割したセクション毎に一定の値を取る関数と定義し、前記第1記憶部内の前記第1時系列データを用いて、前記すべての診断モデルインスタンスの値の二乗を最小化するように前記Bi(Xi)、前記Ci、前記Bj(Xj)、前記Cjを同定することにより、最適化された診断モデルインスタンスをそれぞれ得るモデル最適化部と、
　前記第1の期間と異なる第2の期間において、前記第1～第K制御機器に与えられた前記指示値と、前記第1～第Kセンサによる計測値とを含む第2時系列データを所定時間間隔で記憶する第2記憶部と、
　前記第1～第Kセンサによる最大観測誤差と、前記指示値に対する前記第1～第K制御機器の最大応答誤差とを表す精度情報を記憶する第3記憶部と、
　前記最適化された診断モデルインスタンス毎に、
　前記第1記憶部内の前記第1時系列データを用いて前記最適化された診断モデルインスタンスの値である第1残差を求め、前記第1残差の確率分布である第1確率分布を生成し、
　前記第2記憶部内の前記第2時系列データを用いて前記最適化された診断モデルインスタンスの値である第2残差を求め、前記第2残差の確率分布である第2確率分布を生成し、
　前記最適化された診断モデルインスタンスを生成する元となった2つの制御機器および2つのセンサのそれぞれ毎に、前記第1時系列データにおけるそれぞれに対応する値をそれぞれの前記精度情報に応じて変更し、変更後の時系列データを用いて前記最適化された診断モデルインスタンスの値である第3残差を求め、前記第3残差の確率分布である第3確率分布をそれぞれ生成し、
　前記第2確率分布と前記第1確率分布間の距離を表す第1の拡張KL(Kulback-Leibler divergence)情報量と、前記第2確率分布と各前記第3確率分布のそれぞれとの距離を表す第2の拡張KL情報量とを求め、前記第1および第2拡張KL情報量の合計に対する前記第1拡張KL情報量の比率に応じて判定スコアをそれぞれ計算する判定スコア計算部と、
　前記第1～第K制御機器および前記第1～Kセンサ毎に、前記最適化された診断モデルインスタンスから得られた判定スコアの最小値である総合スコアを求め、第1閾値より大きい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する総合異常判定部と、
　を備える。
　本発明により、個体差のある未知の特性を持つ制御機器の異常検出を高い精度で行うことが可能となる。
本発明の実施の形態としての異常検出装置の構成を概略的に示すブロック図である。 本実施の形態を説明するための事例を説明するための図である。 図2と異なる事例を示す図である。 入力画面の一例を示す図である。 インデックステーブルの一例を示す図である。 データテーブルの一例を示す図である。 目的モデル格納部に格納された複数の目的モデル例を示す図である。 モデル生成部の処理の流れを示すフローチャートである。 モデル生成部の出力であるモデルリストの例を示す図である。 モデル情報の例を示す図である。 モデル最適化部による処理の流れを示すフローチャートである。 図11におけるパラメータ同定演算過程におけるメモリ内容を示す図である。 パラメータ同定演算の様子を示す図である。 パラメータ同定演算の様子を示す図である。 パラメータ同定演算の様子を示す図である。 モデル情報の一例を示す図である。 時系列データ記憶部に記憶された時系列データの例を示す図である。 異常検出部よる処理の流れを示すフローチャートである。 異常スコア計算部による異常スコアの計算処理の流れを示すフローチャートである。 残差および信頼度の計算例を説明する図である。 残差および信頼度の計算例を説明する図である。 異常スコアの計算方法を説明する図である。 精度情報に基づき時系列データを補正する例を示す図である。 総合異常判定部による総合異常判定処理の流れを示すフローチャートである。 計算された異常スコア、および総合スコアの計算例を示す図である。 結果表示の画面例を示す図である。 結果表示の画面例を示す図である。 結果表示の画面例を示す図である。 解析結果画面の例を示す図である。
　図1は本発明の実施の形態としての異常検出装置の構成を概略的に示すブロック図である。
　まず、理解の簡単のため、図1の装置を説明する前に、異常検出を行う対象となる事例について説明する。
　図2は、本実施の形態を説明するための事例を説明するための図である。
　風または空気(媒体)を送り込むファンをダクト内に持たない空調システムを想定する。ダクト本流が1つ存在し、本流ダクトから複数の支流ダクトに分かれる。本例では3つの支流ダクトに分かれる。これら本流および支流は、本発明の媒体経路の一例に相当する。
　各支流ダクトの途中で風量を調節する可変風量装置(VAV：Variable Air Volume)あるいは固定風量装置(CAV: Constant Air Volume)が設置される。
　本例ではVAV1～VAV3が設置される。VAV1～VAV3はそれぞれ異なる箇所に設置される。VAV1～VAV3は、各支流ダクトに流れる風量を調整する弁X1～X3を備える。VAV1～VAV3は、k (k=1..n) 番目の支流の開度指示値x(k)に応じて、弁X1～X3の開きを調整するアクチュエータを含む。弁の開きを調整して支流ダクトに流れる風量を外部に出力する。開度指示値は、図示しないコントローラによって生成され、VAV1～VAV3に与えられる。
　またVAV1～VAV3の近傍に風量センサAF1～AF3が配置される。各支流ダクトの出口が1つの室内に存在し、k (k=1..n) 番目の支流の風量センサは、風量測定値 Q(k)を観測する。風量センサはVAV（制御機器）によって出力制御された風量（出力量）を観測する。
　本装置は、k (k=1..n) 番目の支流のVAVの開度指示値 x(k) と風量測定値 Q(k) を逐次、観測時刻とともに内部の記憶部に記憶し、記憶したデータを元に、アクチュエータ・センサ異常の検出を行う。すなわち、VAV開度指示値x(k)と実際の開度の間に乖離があるか否か(開度異常があるか否か)、風量センサ読み取り値Q(k)と実際の風量の間に乖離があるか否か(風量センサに異常があるか否か)を検出する。
　なおダクト構成は図2に限定されず、様々な態様が考えられる。たとえば図3(A)のように支流ダクトの個数は2個以上の任意の個数nでよい。また図3(B)のように各支流ダクトの出口が複数の室内に存在してもよい。図3(B)の例では、各室内は、解放可能なドアまたは窓によって隣接している。
　以下、図1は異常検出装置について詳細に説明する。
　異常検出装置1は、モデル生成部11、モデル最適化部12、異常検出部13、目的モデル格納部16、基準時系列データ記憶部（第1記憶部）17、モデル格納部18、時系列データ記憶部（第2記憶部）19、結果格納部20、センサ情報記憶部（第3記憶部）21と、条件入力部15、表示部16と、を備える。
　要素11～13は、たとえばプログラムモジュールをCPUに実行させることによって実現される。CPUは、当該プログラムモジュールを、図示しないプログラム記憶部から読み出し、メモリに展開して、実行する。プログラム記憶部は、ハードディスク装置、メモリ装置、CD-R、DVD-R、DVD-RW、フレキシブルディスク等の任意の記憶媒体によって実現できる。
　要素16～21は、ハードディスク装置、メモリ装置、CD-R、DVD-ROM、等の任意の記憶媒体によって構成される。記憶媒体は、異常検出装置1の内部に設けられてもよいし、脱着可能な外部装置でもよい。
　条件入力部15は、ユーザによりデータ入力を行うための入力インターフェースであり、キーボード、マウスなどの任意の入力手段によって構成されることができる。
　表示部16は、ユーザとの出力インターフェースであり、液晶表示装置、プラズマ表示装置、有機EL表示装置、CRT装置などの任意の表示装置によって構成されることができる。表示部16は、たとえばモデル生成部11または異常検出部13から指定されたデータの表示を行う。表示部がタッチパネル機能を有する場合は、条件入力部15は、当該タッチパネル機能を含んでも良い。
　以下、本装置の説明を、モデル生成、モデル最適化、異常検出の順に行う。
[モデル生成]
　図4は、条件入力部15による入力を行うための入力画面の一例を示す。
　入力画面の左側において、建物としてXビル、フロアとして2階～14階のそれぞれの部屋A～Dが選択されている。なお選択するフロア、選択する部屋はそれぞれ1つでもよい。入力画面の右側において、部屋内における異常検出対象となるセンサおよびアクチュエータに対するパラメータ設定を行う。図示の例では、2階～14階のそれぞれの部屋A～Dに対して、共通の設定を行っている。ここでは図2の事例を想定し、3つの支流ダクトが存在し、各支流ダクトに配置されたVAV1～VAV3に対し、VAV開度指示値と、風量をそれぞれ観測パラメータとして選択している。
　条件入力部15による以上の入力により、以下の形式を有するポイントオブジェクト名が設定される。
<ビルID.フロアID.ポイント大種別.ポイント小種別.ポイントID.状態量>
　ポイントオブジェクト名は、区切り文字(ここでは “.” )によって区分された6つのフィールドを含む。
　ビルIDは、選択した建物のIDである。　
　フロアIDは、選択した階および部屋を特定するIDである。　
　ポイント大種別は、部屋内での場所種別である。たとえば「I」(Interior)は部屋の内部、「P」(Perimeter)は窓際を表す。　
　ポイント小種別は、機器種別である。たとえばVAV、CAVがある。　
　ポイントIDは、ポイント小種別で指定された機器のIDである。　
　状態量は、当該指定された機器に対して観測される物理量または指示値を示す。
　図4の例では、たとえばXビルの12階の部屋Cに対して、1番目の支流ダクト(あるいは1番目のVAV)に対して、
<Xビル.12C.I.VAV.1.開度指示値>と、<Xビル.12C.I.VAV.1.風量>が設定される。　
　また、2番目の支流ダクト(あるいは2番目のVAV)に対して、
<Xビル.12C.I.VAV.2.開度指示値>と、<Xビル.12C.I.VAV.2.風量>が設定され、
　3目の支流ダクト(あるいは3番目のVAV)に対して、
<Xビル.12C.I.VAV.3.開度指示値>と、<Xビル.12C.I.VAV.3.風量>が設定される。
　図1の基準時系列データ記憶部（第1記憶部）17は、インデックステーブルと、データテーブルをあらかじめ記憶している。
　図5はインデックステーブル171の一例を示す。　
　インデックステーブル171は、ハンドルIDと、ポイントオブジェクト名から構成される。インデックステーブル171は、たとえば図4の入力画面で設定可能なすべてのポイントオブジェクト名とそのハンドルIDとを記憶している。
　図6はデータテーブル172の一例を示す。　
　データテーブル172は、インデックステーブル171のハンドルIDごとに、該当するポイントオブジェクト名で特定されるセンサ・アクチュエータの観測値の履歴を、基準時系列データ（第1時系列データ）として記憶している。データテーブル171に記憶された基準時系列データ（第1時系列データ）は、いずれもセンサ・アクチュエータが正常であるときの第1の期間に取得されたものである。各ハンドルIDに対応する観測値は、それぞれ一定時間間隔(本例では10分ごと)に、同時に記録されている。
　目的モデル格納部16には、ビル空調システムにおいて下記の数式1が適用可能な目的モデルを、複数格納している。目的モデルは診断対象システムにおける部分構造(本例ではビル空調システムにおける任意の2つの制御機器間（あるいは支流ダクト間)）に適用可能なものである。
　数式1に含まれる (Ci + βi(Xi) ) × (Qi/Qj)^2－ (Cj + βj(Xj) ) が、本発明の目的モデルに対応する。なお「^」はべき乗を意味する。数式1は、目的モデルのCi, Cj, βi, βjを、当該目的モデルの値(e)の二乗(R) を最小にするように求めることを示している。　
(数式 1)
min_{Ci, Cj, βi, βj}　
R　=　{ (Ci + βi(Xi) ) × (Qi/Qj) ^2 － (Cj + βj(Xj) ) }^2, ただしCi,Cj>0
Ci, Cj ：定数項、βi, βj ：未知特性項、Xi：第1の観測変数、Qi：第2の観測変数,Xj：第3の観測変数,Qj：第4の観測変数
　なお数式1をより一般化した以下の式1aを用いてもよい。式1aのf(Y)は第2および第4観測変数の演算項であり、既知の形の関数である。f(Y)はf(Qi,Qj)と記述されてもよい。この場合、(Ci + βi(Xi) ] ) × f(Y) － (Cj + βj(Xj) )が目的モデルに対応する。　
(数式 1a)
min_{Ci, Cj, βi, βj}　
R　=　{ (Ci + βi(Xi) ] ) × f(Y) － (Cj + βj(Xj) ) }^2
　上記の数式1に含まれる目的モデルは任意の2つの支流ダクト内に流れる媒体(流体または気体)のエネルギー保存則をモデル化したものである。このエネルギー保存則は以下の数式1bに示される。図2のような3つのVAVが配置されたダクト構成を例に、目的モデルの導出過程を簡単に説明する。
Figure JPOXMLDOC01-appb-M000001
　Ki1,Ki2,Kj1,Kj2は、対象ゾーンのインテリイア/ペリメータエリアに設置されたダクトの圧力損失係数を表す定数である。
　i,jは異なる2つの支流を示し、ダクト(またはVAV)のIDである。
　Ki1,Kj1は、各支流への分岐点から各支流端に至る管路で定常的に損失する圧力損失係数(支流ダクトの固有の圧力損失係数)である。
　Ki2,Kj2は各支流ダクトのVAVの開度に応じて可変的に損失する各支流ダクトの圧力損失係数を示す。
　Qi,Qjは風量センサが示す風量計測値(m3/s)である。
　vi,vjはダクト風量を調整する弁の開度指示値(0～1の実数で単位無し)である。
　数式1bは、支流iのVAV指示値(VAV開度)とVAV風量、ならびに支流jのVAV指示値(VAV開度)から、支流jのVAV風量を推定できことを意味している。この関係式はVAVの数が変更しても成立するとともに、外乱の影響を受けにくいため汎用性が高い。
　これらのパラメータを用いて、ベルヌーイの定理から、任意の2つの支流ダクトに対して、数式1bが導出される。ベルヌーイの定理に基づく数式1bの導出過程は、説明の簡略化のため、割愛する。この数式1bを式変形することで、上述の一般化された目的モデルを導出できる。
Figure JPOXMLDOC01-appb-M000002
　図7は、目的モデル格納部16に格納された複数の目的モデル例を示す。　
　F-01,F02,F03,F04,F05,F06….の目的モデル(エントリ)が格納されている、各目的モデルに含まれる< >内に記述される文字列は、ポイントオブジェクト名に対する抽象表現を意味している。
　[s1]、 [s2] 、 [s3] 、 [s4] および [s5]の変数 はそれぞれ任意の文字列に対応づけられる。ただし、1つのエントリに同じ変数が複数存在する場合、各変数は同じ文字列に対応づけられなければならない。たとえばF-01の目的モデルには、［s1］が3つ含まれるが、これらの［s1］はそれぞれ同一の文字列に対応づけられなければならない。
　モデル生成部11は、入力情報解析部111とインスタンス生成部112とを備える。モデル生成部11は、これらの要素111，112を用いて、条件入力部15で入力された情報(ポイントオブジェクト名)と、目的モデル格納部16に基づき、数式1のインスタンス生成を行う。
　ここにおけるインスタンスとは、数式1の任意の観測変数Xi, Xj, Qi,Qj に対して、具体的な観測変数を割り当てることを意味する。たとえば数式1のXiに “Xビル.2A.I.VAV.1.開度指示値” 、Qiに “Xビル.2A.I.VAV.1風量” を割り当てるなどである。
　以下、図8を参照して、モデル生成部11の処理の詳細を説明する。
　図8は、モデル生成部11の処理の流れを示すフローチャートである。　
　モデル生成部11は、条件入力部15から利用者が情報を入力したことをトリガーとして処理を開始する。
　入力情報読込部111が、条件入力部15により登録された情報(複数のポイントオブジェクト名)と、基準時系列データ記憶部17のインデックステーブル171をメモリに読み込む(S11)。
　各ポイントオブジェクト名における上位3つの情報(ビルID、フロアID、ポイント大分類ID)を抽出し、抽出した上位3つの情報を含む要素のリスト(部分構成リスト)を作成する(S12)。
　図4の入力事例の場合、{“Xビル.2A.I”, “Xビル.2B.I”, …, “Xビル.14C.I”, “Xビル.14D.I”} のような52個の要素を含む部分構成リストが作成される。
　次に、部分構成リストにおける各要素に対して、以下のステップS13～S16を実行する。
　ステップS13で、部分構成リストから、まだ選択されていない要素を選択する。
　ステップS14では、ステップS13で選択された要素に対して、診断モデルインスタンスを生成する(S14)。
　たとえば選択された要素 が“Xビル.2A.I”であるとすると以下のようになる。
　条件入力部15から入力されたポイントオブジェクト名のうち、この要素“Xビル.2A.I”を含むポイントオブジェクト名を列挙すると、以下の6個のポイントオブジェクト名が得られる。　
　　　Xビル.2A.I.VAV.1.開度指示値
　　　Xビル.2A.I.VAV.1.風量
　　　Xビル.2A.I.VAV.2開度指示値
　　　Xビル.2A.I.VAV.2.風量
　　　Xビル.2A.I.VAV.3.開度指示値
　　　Xビル.2A.I.VAV.3.風量
　これらのポイントオブジェクト名を、目的モデル格納部16に登録されている各モデルと照合すると、下記のF-01のモデルにマッチすることが分かる。照合は文字列マッチングで行えばよい。つまり、[s1],[s2]等の各変数に矛盾が発生しない(同じ変数には同じ文字列が入る)ようなモデルを選択する。　
　F-01: (C1 + β1(<[s1].[s2].[s3].VAV.[s4].開度指示値>)) 
　　　　*(<[s1].[s2].[s3].VAV.[s4].風量>/<[s1].[s2].[s3].VAV.[s5].風量>) 2 -(C2 + β2(<[s1].[s2].[s3].VAV.[s5].開度指示値>)) 
　このF-01のモデルにおいて、[s1]～[s3]は一意に決まるが、[s4], [s5]に関しては6通りの組合せが発生する。
　具体的に、[s1]は“Xビル”、[s2]は“2A”、[s3]は“I” で一意に決まるが、[s4], [s5]の組み合わせは、以下の6通り存在する。なお、［s4］,[s5]はいずれもポイントIDに対応する。　
　　　([s4], [s5]) =(1, 2), (2, 1), (1, 3), (3, 1), (2, 3), (3, 2)
　よって、F-01のモデルに対して、下記の6通りのインスタンス(診断モデルインスタンス)<1>～<6>が得られる(図9参照)。
Figure JPOXMLDOC01-appb-T000003
　ステップS15では、ステップS14で生成した各診断モデルインスタンスに対応して、目的モデルに照合するための部分診断リスト(パラメータテーブルとアクセステーブル)を登録する。
　ここでパラメータテーブルは、各診断モデルインスタンスの係数(Ci, Cj, βi, βj)を格納したテーブルである。パラメータテーブルの例を図9の右に示す(「PARMS」属性にリンクされている)。
　たとえば上記表の最下のインスタンス<6>に対しては、F-01の係数C1として係数C［3］、F-01の係数C2として係数C［2］、F-01の係数β1として係数β［3］、F-01の係数β2としてβ[2]が格納される。すなわち、当該最下のインスタンスに対して、F-01のC1,C2,β1,β2(あるいは数式1内のCi,Cj,β2, β2)に対応して、“C[3] ,C[2],β[3],β[2]”がパラメータテーブルに登録される。なお現時点ではこれらのパラメータの値は同定されていない(未知である)。
　同様にして他の残りの5つのインスタンス<1>～<5>に対しても同様にパラメータが登録される。
　アクセステーブルは、各診断モデルインスタンスのそれぞれに含まれる観測変数(開度指示値、風量)に対応するIDを格納する。アクセステーブルの例が図9の右に示される(「ACCESS」属性にリンクされている)。
　アクセステーブルは、基準時系列データ記憶部17に記憶された基準時系列データへのマッピング、および後述する異常診断において用いる時系列データ記憶部19に記憶された時系列データへのマッピングのために用いられる。　
　たとえばインスタンス<6>の場合、
　<Xビル.12C.I.VAV.3.開度指示値>
　<Xビル.12C.I.VAV.2.開度指示値>
　<Xビル.12C.I.VAV.3.風量>
　<Xビル.12C.I.VAV.3.風量>
　の4つの観測変数が存在し、
　それぞれ対応するハンドルIDは、135,134,145,144である(図5のインデックステーブル171を参照)。したがって、上記インスタンス<6>に対して、“135,134,145,144”がアクセステーブルに登録される。
　以上に説明したステップS14,S15の処理を、部分構成リストの残りのすべての要素についても同様に行う。
　部分構成リストに含まれるすべての要素についての処理が終了すると、ビルの階およびフロア毎のモデル情報を含むモデルリストが、モデル格納部18に格納される(ステップS17)。
　モデルリストはモデル生成部の出力であり、図9に示される。
　図9の右には、「12C」(12階の部屋C)に対して生成されたモデル情報の例が取り出して示されている。他の残りの階および部屋についても同様の形式のモデル情報が存在する。これらのモデル情報の集合がモデルリストに相当する。
　図9に示すように、モデルリストは、条件入力部で選択された各階および各部屋のモデル情報を含み、1～Nmの番号がモデル情報の番号である。ビルIDとフロアIDの組ごとに番号が設定され、モデル情報が生成されるごとに順次1～Nmの番号が増える。図10に「2A」(2階の部屋A)に対して生成されたモデル情報の例を示す。「2A」のモデル情報は上記の処理フローで最初に作成されるため番号1が付与され、この時点ではまだ他のモデル情報は存在しない。「12C」に対するモデル情報の生成は39番目に行われ、番号39が付与されている(図9参照)。なお、本例では最後の番号Nmは52となる。
　ここでモデル情報はそれぞれ以下の属性を含む。
“ID”、“Fig”、“BIND”、“PARMS”、“ACCESS”、“C”、“segment”,“B”,“w”
　たとえば、Xビルの12階の部屋Cのモデル情報の番号は39番であり、“ID”は、モデル情報の識別子“Xビル.12C”へのリンクを示す。
　“Fig”は、条件入力部15による入力画面の画面データが保存されたファイルへのリンクを示す。
　“Bind”は、前述した診断モデルインスタンスへのリンクを示す。
　“PARMS”は、先に少し説明したが、パラメータテーブルにリンクされ、“ACCESS”がアクセステーブルにリンクされる。
　属性“C”(定数項),“segment”(セグメント)、“B”(未知特性項)、“w”(信頼値)に対しては現時点ではまだ何にも設定されていない。詳細は後述する。
　なお図4の入力画面ではポイント大種別が“I”の場合を示したが、1つの部屋に“I”と“P”が混在する場合もある。この場合は、“I”と“P”それぞれごとに別個に処理が行われる(上述したように目的モデルへのマッチングは、ポイントオブジェクト名における上位3つの情報(ビルID、フロアID、ポイント大分類ID)により行われる)。したがって、この場合、1つの部屋(たとえば12階の部屋C)に対して、モデル情報が、“I”と“P”のそれぞれごとに生成される。
[モデル最適化]
　図1のモデル最適化部12は、データ読込部121、未知特性項初期演算部122、定数項収束演算部123、未知特性項収束演算部124を備える。
　モデル最適化部12は、モデル情報毎に、基準時系列データ記憶部17に記憶された基準時系列データに対して、診断モデルインスタンスの値の二乗(数式1のR)を最小化するように当該診断モデルインスタンスのパラメータを学習する。
　データ読込部121は、基準時系列データ記憶部17に記憶された基準時系列データをメモリに読み込む。
　未知特性項初期演算部122は、各診断モデルインスタンスに含まれる定数項を固定するとともに未知特性項を共通とする条件下で、メモリに読み込んだ基準時系列データに対して誤差(診断モデルインスタンスの値)の二乗が最小になるように、未知特性項を同定し、メモリに記憶する。
　定数項収束演算部123では、現在メモリに記憶されている状態で未知特性項を固定した上、メモリに読み込んだ基準時系列データに対して誤差の二乗が最小になるように定数項を同定しメモリに記憶する。
　未知特性項収束演算部124は、現在メモリに記憶されている状態で定数項を固定した上、メモリに読み込んだ基準時系列データに対して誤差の二乗が最小になるように、未知特性項を同定し、メモリに記憶する。
　定数項収束演算部123と未知特性項収束演算部124は、終了条件を満たすまで交互に処理を繰り返し実行し、終了条件を満たした段階で、そのときのパラメータをモデル格納部18に格納する。このように求められたパラメータを設定した診断モデルインスタンスは、最適化された診断モデルインスタンスと称される。
　図11は、モデル最適化部12による処理の流れを示すフローチャートである。モデル最適化部12は、モデル生成部11の後に実行される。したがって、モデル最適化部12の実行前には、Xビルのユーザ選択された階・部屋に対するモデル情報がモデル格納部18に記憶されている。
　まず、データ読込部21により以下のステップS21～S23の処理を実行する。
　モデル格納部18に記憶されたモデルリストにおいて任意の番号のモデル情報をメモリに読み込む (S21)。たとえば39番の「Xビル、12階、部屋C」のモデル情報を読み込む。
　読み込んだモデル情報に含まれる診断モデルインスタンス毎に、ACCESS属性を利用して、基準時系列データ記憶部17の基準時系列データから該当するデータ(ACCESS属性のハンドルIDに対応する観測値と、タイムスタンプ(時刻情報))をメモリに読み込む(S22)。
　ステップS22で読み込んだデータに対して欠損値処理を行う(S23)。欠損値処理は、既存の方法で行う。たとえば、入力されたデータに対して大きな値変更をすることがなく、欠損値や外れ値の除去能力が比較的高く、かつ計算量が少ない修正メジアンフィルタを適用することが望ましい。
　次に、未知特性項初期演算部122、定数項収束演算部123および未知特性項収束演算部124によりステップS24～S28を行う。
　ステップS24は未知特性項初期演算部122により行う。まず、未知特性項であるβk(X) は既定の関数で与えられていないため、βk(X)を、α[k, s(X)] と置き換える。α[k, s(X)]は、Xの取り得る値の範囲を一定幅で分割したM個のセグメントのそれぞれ毎に定数を与える関数である。
　すなわち、s(X)は、Xが与えられた場合に、Xが属するセグメントのインデックスを返す関数であり、与えられたXに対して1～Mのいずれかの自然数を出力する(M：セグメント数)。kは、数式1の例ではiまたはjに相当する未知特性項識別子である。したがって、α[k, s(X)]は、Xが与えられたとき、Xが属するセグメントの定数を返す。
　上記の置き換えを用いることにより、最適化を示す式である数式1は、次の数式2で表すことができる。　
(数式 2)
min_{C[i], C[j], α[i,1], …,α[i,M], α[j,1], …, α[j,M]}　
 R　=　{ (C[i] + α[i, s(Xi) ] ) × (Qi/Qj)^2 － (C[j] + α[j, s(Xj) ] ) }^2,　i≠j
　この数式2は Rを最小化するようにC[i], C[j], α[i,1], …,α[i,M], α[j,1], …, α[j,M] の計(2×M + 2)個のパラメータを求める(最適化する)ことを意味する。C[i]>0, C[j]>0を制約条件とする。
　この問題を一度に求解するのは困難である。このような問題を解く場合の1つとして、固定する目的変数と、最適化する目的変数とを交互に入れ替えながら反復演算で求解する方法が知られている。ここでは、以下の数式3、数式4、数式5による3つのステップの反復演算によって、数式2を求解する。　
(数式 3)
min_{α[0,1], …, α[0,M] } 
R　=　{ (C0 + α[0, s(Xi) ] ) * (Qi/Qj)^2 － (C0 + α[0, s(Xj) ] ) }^2

(数式4)
min_{C[i], C[j]}　
　　　R　=　{ (C[i] + α[i, s(Xi) ] ) * (Qi/Qj)^2 － (C[j] + α[j, s(Xj)] ) }^2

(数式5)
min_{α[i,1], …,α[i, M], α[j, 1], …, α[j, M]} 
　　 R　=　{ (C[i] + α[i, s(Xi) ] ) * (Qi/Qj)^2 － (C[j] + α[j, s(Xj) ] ) }^2

(ステップ1)1つめのステップは、未知特性項初期演算部122において、数式3のインスタンスである下記の数式3.1を求解する(S25)。
(数式3.1)
min_{α[0,1], …, α[0,M]}　
R　=　{ (C0 + α[0, s(xi) ] ) × (Qi/Qj)^2 － (C0 + α[0, s(xj) ] ) }^2,　i≠j
　以下、事例に従って (数式3.1)を説明する。
　3つの各支流の固有のダクト圧力損失係数に相当するC1, C2, C3 が同じ値C0( = 1)を持つと仮定する。すなわち各支流のダクトの長さ・形状・材質が同じと仮定する。なおC1, C2, C3は、VAV１～３が配置された箇所（ダクト）に応じた未知の定数項であるといえる。
　さらに、VAV1～VAV3のそれぞれの開度指示値から推定されるダクト圧力損失係数に相当するα[1,1:M], α[2,1:M], α[3,1:M] も、機器ごとの個体差がなく共通の支流ダクト圧力損失係数α[0,1:M] (すなわちα[0,1], …, α[0,M])を持つと仮定する。つまり、α[1,1], α[2,1], α[3,1]はいずれも同じ値α[0,1]をもち、α[1,2], α[2,2], α[3,2]はいずれも同じ値α[0,2]をもち、・・・α[1,M], α[2,M], α[3,M]はいずれも同じ値α[0,M]をもつと仮定する。
　このとき、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、共通の支流ダクト圧力損失係数に相当するα[0,1:M]を求める。
　
　具体的には、基準時系列データ記憶部17から読み込んだデータを、 診断モデルインスタンス<1>～<6>毎に、数式3.1に代入することで、下記のように、データ数(=n)×6の非線形連立方程式(数式3.2)が立てられる。これを既存のソルバで解くことで、α[0,1:M] (=α[0,1], …,α[0,M])が求まる(S24～S25)。すなわち各式の左辺が同時にできるだけゼロに近づくように最適化することで、α[0,1:M] (=α[0,1], …,α[0,M])が求まる。数式3.2において、R[a,b]の“a”は、診断モデルインスタンスを識別する番号であり、“b”はn個の各データを識別する番号である。　
 (数式3.2)
{ (C0 + α[0, s(x1[1]) ] ) × (Q1[1] / Q2[1] )^2 － (C0 + α[0, s(x2[1]) ] ) }^2 －R[1,1] = 0
…
{ (C0 + α[0, s(x1[n]) ] ) × (Q1[n] / Q2[n] )^2 － (C0 + α[0, s(x2[n]) ] ) }^2 － R[1,n] = 0
{ (C0 + α[0, s(x1[1]) ] ) × (Q1[1] / Q3[1] )^2 － (C0 + α[0, s(x3[1]) ] ) }^2 － R[2,1] = 0
…
{ (C0 + α[0, s(x1[n]) ] ) × (Q1[n] / Q3[n] )^2 － (C0 + α[0, s(x3[n]) ] ) }^2 － R[2,n] = 0
{ (C0 + α[0, s(x2[1]) ] ) × (Q2[1] / Q1[1] )^2 － (C0 + α[0, s(x1[1]) ] ) }^2 － R[3,1] = 0
…
{ (C0 + α[0, s(x2[n]) ] ) × (Q2[n] / Q1[n] )^2 － (C0 + α[0, s(x1[n]) ] ) }^2 － R[3,n] = 0
・　　・・・
{ (C0 + α[0, s(x3[1]) ] ) × (Q3[1] / Q2[1] )^2 － (C0 + α[0, s(x2[1]) ] ) }^2 － R[6,1] = 0
…
{ (C0 + α[0, s(x3[n]) ] ) × (Q3[n] / Q2[n] )^2 － (C0 + α[0, s(x2[n]) ] ) }^2 － R[6,n] = 0
　図13に、本ステップ1(S25)の演算の様子を示す。　
　図13の左に示すように、C1, C2, C3 が同じ値C0( = 1)として固定され、またα[0,1:M] (=α[0,1], …,α[0,M])の初期値を1としている。そして、上記数式3.2を解くことで、図13の右のように、α[0,1:M] (=α[0,1], …,α[0,M])が求められている。X軸がVAV開度指示値で、Y軸がα[0,s(xk)]+C0である。
　図12は、ステップ1～3(S25～S27)の演算過程におけるメモリ内容を示すもので、最適化する変数はメモリ上で配列テーブルとして保持されている。本ステップ1のメモリ内容は図12の上に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ1では上述のようにC1, C2, C3 が同じ値C0( = 1)を持つとして固定される(白抜きセル)。α[1,1:M], α[2,1:M], α[3,1:M]は共通のα[0,1:M]を持つと仮定した上で、α[0,1:M]が最適化される(網掛けのセル)。
(ステップ2)2つめのステップは、定数項収束演算部123において、数式4のインスタンスである下記の数式4.1を求解する(S26)。　
(数式4.1)
min_{C[i], C[j]}　
　　　R　=　{ (Ci + α[i, s(xi) ] ) ×(Qi/Qj)^2 － (Cj + α[j, s(xj) ] ) }^2, i≠j
　以下、事例に従って(数式4.1)を説明する。　
　VAV1～VAV3の開度に対する圧力損失係数に相当するα[1,1:M], α[2,1:M], α[3,1:M]を現在の状態で固定(前述のように最初は、α[1,1:M], α[2,1:M], α[3,1:M]はすべて共通のα[0,1:M])し、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、各支流のダクト圧力損失係数に相当する C1, C2, C3を求める。
　ステップ1と同様に、診断モデルインスタンス<1>～<6>毎に、読み込んだデータを入力することで、データ数(=nとする)×6の連立方程式が立てられ、これを既存のソルバで解くことでC1,C2,C3が求まる (S26)。
　図14に、本ステップ2(S26)の演算の様子を示す。　
　図左のC1old,C2 old,C3 oldは、本ステップによる最適化前のC1,C2,C3であり、図右のC1new,C2 new,C3 newは、本ステップによる最適化後のC1,C2,C3である。
　各支流ダクトに対応するα[k,s(xk)](すなわちα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)])は、本ステップでは固定される。
　なお、X軸は、VAV開度指示値であり、具体的にx1はVAV1の開度指示値、x2はVAV2の開度指示値、x3はVAV3の開度指示値である。Y軸は、α[k,s(x)]+Ckである。
　本ステップ2(S26)でのメモリ内容は図12の上から2番目に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ2では上述のようにα[1,1:M], α[2,1:M], α[3,1:M]を固定しつつ、C1,C2,C3を最適化している。
(ステップ3)3つめのステップは、未知特性項収束演算部124において、数式5のインスタンスである下記の数式5.1を求解する。　
(数式5.1)
min_{α[i,1], …,α[i,M], α[j,1], …, α[j,M]}　
R　=　{ (C[i] + α[i, s(xi) ] ) ×(Qi/Qj)^2 － (C[j] + α[j, s(xj) ] ) }^2,　 i, j=1,..,3, i≠j
　以下、事例に従って(数式5.1)を説明する。　
　各支流のダクト圧力損失係数に相当する C1, C2, C3を現在の状態(S26で求めた値)で固定し、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、VAV1～VAV3の開度指示値に対するダクト圧力損失係数に相当するα[1,1:M] (= α[1,1], …,α[1,M])、α[2,1:M](=α[2,1], …, α[2,M])、α[3,1:M](=α[3,1], …, α[3,M])を求める。具体的には、これまで同様に、診断モデルインスタンス<1>～<6>毎に、読み込んデータを入力することで、データ数(=nとする)×6の連立方程式を立て、これを既存のソルバで解くことで、α[1,1:M]、α[2,1:M]、α[3,1:M]が求まる (S27)。
　図15に、本ステップ3(S27)の演算の様子を示す。
　図左のαold [1,s(x1)]、αold [2,s(x2)]、αold [3,s(x3)]は、本ステップによる最適化前のα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)]であり、αnew [1,s(x1)]、αnew [2,s(x2)]、αnew [3,s(x3)]は本ステップによる最適化後のα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)]である。
　各支流ダクトに対応するC1,C2,C3は、本ステップでは固定される。
　なお、図右の網掛けのブロックは最適化の前後でα[k,s(x)]の値が変化したことを示し、網掛けされていないブロックは最適化の前後でα[k,s(x)]の値が変化しなかったことを示す。
　本ステップ3(S27)でのメモリ内容は図12の下から2番目に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ3では上述のようにC1,C2,C3を固定しつつ、α[1,1:M], α[2,1:M], α[3,1:M]を最適化している。
　ここでステップS27では後述の異常判定の信頼性向上のため以下の処理を行っても良い。
　ステップS22でメモリに読み込んだデータに含まれるVAV開度指示値の同時分布に疎な部分がある場合、疎な部分のデータにより求められる部分解 α[k,s] は、真に求めるべき値に対して誤差が大きくなる可能性がある。誤差が大きくなると、後述の異常検出部13における異常判定に誤りが生じやすくなる。そこで、後の異常判定においてαの信頼性に基づいた処理を行うことができるように、α[k,s]に対する信頼値 w[k,s]を設定する。
　具体的に、α[1,1:M]( =α[1,1]、α[1,2] , …,α[1,M])に対し、信頼値w[1,1:M](= w [1,1]、w [1,2] , …, w [1,M]))を計算し、α[2,1:M]( =α[2,1]、α[2,2] , …,α[2,M])に対し、信頼値w[2,1:M](= w [2,1]、w [2,2] , …, w [2,M]))を計算し、α[3,1:M]( =α[3,1]、α[3,2] , …,α[3,M])に対し、信頼値w[3,1:M](= w [3,1]、w [3,2] , …, w [3,M]))を計算する。
　この際、区間s(=1,2,…,M)に含まれるデータ数が多いほど、信頼値が高くなるようにする。たとえばデータ数が多いほど高い信頼値を与える関数またはテーブルを用意しておき、当該関数またはテーブルに、該当区間に含まれるデータ数を与えることにより、当該区間の信頼値を計算する。
　本ステップS27で信頼値の計算を行う場合のメモリ内容が図12の下に示される。なお信頼値の計算は、一度行えばよいため、ステップS27の実行のたびに、繰り返し行う必要はない。
　以上に説明したステップ2，3(図11のS26～S27)を、終了条件を満たすまで繰り返し実行することで(S28)、環境に対応した、固有の定数項Ck (k=1,2,3)(例：支流ダクトの固有の圧力損失係数)と、機器特有の未知特性項α[k, s(xk)](k=1,2,3, s(xk)=1,2,3,…M) (例：VAV開度指示値から推定される圧力損失係数)を求めることができる。
　ここで終了条件は、利用したソルバに応じて決定されるものであるが、たとえば、「解に収束」「目的変数の変化量が許容範囲より小さい」「繰り返し計算回数が最大数に到達」が挙げられる。
　終了条件が満足されたら、求められたモデルパラメータ情報として、Ck (k=1,2,3)、α[k, s(xk)](k=1,2,3, s(xk)=1,2,3,…M))を、モデル格納部18に記憶させる(S29)。
　具体的に、モデル格納部18のモデル情報の属性“C”(定数項),“segment”(セグメント)、“B”(未知特性項)、“w”(信頼値)にリンクさせて、Ck、s(xk), α[k, s(xk)],w[k, s(xk)]をそれぞれ記憶させる。
　ステップS29でパラメータ情報が追加された後のモデル情報の一例を図16に示す。
[異常検出]
　時系列データ記憶部(第2時系列データ記憶部)19は、異常検出の対象とする機器の状態量(指示値、風量など)を、一定時間間隔で順次、時系列データ(第2時系列データ)としてリルタイムに記憶する。具体的に、図5に示したような各ハンドルIDに対応する状態量を逐次記憶する。時系列データ記憶部19に記憶された時系列データの例を図17に示す。時系列データ記憶部19がデータを記憶する期間は、基準時系列データ記憶部17が記憶する第1の期間と異なる第２の期間である。第２の期間は第1の期間より後でもよいし、前でもよい。
　センサ情報記憶部21は、センサ・アクチュエータの精度情報を記憶する。精度情報は、たとえば工場出荷でメーカが保証している精度を用いる。本例では風量センサ毎の最大検出誤差が、たとえば±5％、または±0.5(m/s)などとして、センサ情報記憶部21に記憶させている。また開度指示値に対するアクチュエータの最大応答誤差がセンサ情報記憶部21に記憶されている。アクチュエータの最大応答誤差は、指示値に対して最大でどの程度、弁の実動作開度に誤差が生じるかを示す。なおアクチュエータの精度情報は、開度指示値での表現に換算して、センサ情報記憶部21に記憶されてもよい。たとえば弁度の開度誤差範囲に対応する開度指示値範囲を求め、求めた範囲を精度情報として記憶してもよい。
　異常検出部13は、モデル読込部131、データ選択部132、異常スコア計算部(判定スコア計算部)133、総合異常判定部134、異常検出提示情報生成部135を備える。
　異常検出部13は、時系列データ記憶部19から診断対象とする任意の期間分の時系列データを読み込み、読み込んだ時系列データから制御機器(アクチュエータ)・センサに異常があるか否かを判定する。そして、判定結果を結果格納部20へ格納するとともに、表示部16へ表示する。
　これまで説明したモデル生成部11、モデル最適化部12が、異常検出装置1の導入直後や、異常検出の対象となる機器を取り巻く環境に変化があった場合に実行されるブロックであるのに対し、異常検出部13は、常時、もしくは、定期間隔(例えば毎日1回など)で実行される。異常検出部13の処理の実行前には、モデル生成部11およびモデル最適化部12によって生成および更新されたモデル情報がモデル格納部18に格納されていることが前提となる。
　以下では、1日1回、VAV開度とVAV風量の異常診断を実行する場合を例にして説明を行う。
　図18は、異常検出部13における要素131～134による処理の流れを示すフローチャートである。
　まず、モデル読込部131がモデル格納部18に記憶されたモデル情報を読み込む(S301)。より詳細に、ビルIDをキーにして、モデル格納部18に格納されているモデル情報(各階の各部屋)を読み込む。モデル情報毎に、以降のステップを行って、異常診断を行う。説明の簡単のため、12階の部屋Cのモデル情報を想定して説明を行う。
　データ選択部132が、モデル情報の診断モデルインスタンス毎に、ACCESS属性を利用して、時系列データ記憶部19から、対応する時系列データ(ACCESS属性に対応する値と時刻)をメモリに読み込む(S303)。
　データ選択部132は、読み込んだ時系列データに対して欠損値処理を行う(S304)。欠損値処理は既存の方法を用いればよい。読み込んだ時系列データに対して大きな値変更をすることがなく、欠損値や外れ値の除去能力が比較的高く、かつ計算量が少ない修正メジアンフィルタを適用することが望ましい。本ステップは省略することも可能である。
　次に、異常スコア計算部133による異常スコア(判定スコア)の計算を行う(S305)。この処理を、当該モデル情報における診断モデルインスタンス毎に繰り返し実行する。
　図19は異常スコア計算部133による異常スコアの計算処理の流れを示すフローチャートである。
　まず、最初の診断モデルインスタンスに対応するPRMS属性、ACCESS属性に基づき、該当するパラメータ情報と、時系列データの時刻毎のデータとを参照し、数式6(残差の計算式)、および数式7(信頼度の計算式)を計算することで、残差ベクトルと信頼度ベクトルを求める(S305a,S305b)。以下図20および図21を用いて、詳細を説明する。　
(数式6)
　　　 e =　( C[i] + α[i, s(xi) ] ) ×( Qi / Qj )^2 － ( C[j] + α[j, s(xj) ] ) 
(数式7)
　　　 r= w[i,s(xi)] × w[j, s(xj)]
　
　たとえば1番目の診断モデルインスタンス<1>の場合、まず時系列データの先頭データに対して、図20に示す参照を行う。2010/8/22/00:10の時系列データ(診断期間の一番始めの時系列データ)のうち、ACCESS属性に示されるID133,134,143,144の値が参照される。ID133はVAV1の開度指示値、ID 134はVAV2の開度指示値、ID 143はVAV1の風量、ID 144はVAV2の風量である(図5参照)。
　またVAV1、VAV2の定数項は、C1=2.31、C2=1.99である。
　また、VAV1の未知特性項B1(α[1, s(x1) ])、VAV2の未知特性項B2(α[2, s(x2) ])の値は、それぞれ、0.43、0.44である。すなわち、VAV1の観測変数値(開度指示値)は、0.71(ID133の時系列データ値)であり、当該0.71はsegment属性のセグメント「0．8」に属するから、α[1, s(0.71) ]=α[1, 0.8 ]=0.43となる。同様にVAV2の観測変数値(開度指示値)のは、0.71(ID134の時系列データ値)であり、当該0.71はsegment属性のセグメント「0．8」に属するから、α[2, s(0.71) ]=α[2, 0.8 ]=0.44となる。
　また各未知特性項B1,B2の信頼値は、segment属性のセグメント「0．8」に対応する0.98である。
　また風量センサAF1,AF2の風量値Q1,Q2は、0.20,0.21(それぞれID143,144の時系列データ値)である。
　したがって数式6により、診断モデルインスタンス<1>の残差は以下のように計算される。　
　　　 e =　( C[i] + α[i, s(xi) ] ) ×( Qi / Qj )^2 － ( C[j] + α[j, s(xj) ] ) 
　　　　= ( C[1] + α[1, s(x1) ] ) ×( Q1 / Q2 )^2 － ( C[2] + α[2, s(x2) ] ) 
　　　　= (2.31 + α[1, s(0.71) ] ) ×( 0.20 / 0.21 )^2 － ( 1.99 + α[2, s(0.71) ] ) 
　　　　 = (2.31 + α[1, 0.8 ] ) ×( 0.20 / 0.21 )^2 － ( 1.99 + α[2, 0.8 ] )
　　　　= (2.31 + 0.43 ) ×( 0.20 / 0.21 )^2 － ( 1.99 + 0.44 ) 
　　　　 =0.06
　また数式7により、診断モデルインスタンス<1>の信頼度は以下のように計算される。
　　　r= w[i,s(xi)] × w[j, s(xj)]
　　　 = w[1,s(x1)] × w[2, s(x2)]
　　　　= w[1,s(0.71)] × w[2, s(0.71)]
　　　 = w[1,0.8] × w[2, 08]
　　　　= 0.98 × 0.98
　　　 =0.96
　以上に示した数式6および数式7の計算を、時系列データの2番目以降の時刻のデータについても同様にして行う。2番目の時刻のデータの場合の参照を示すと、図21のようになる。この場合、数式6によりe=-0.03、数式7によりr=0.88が計算される。
　以上の処理により、診断モデルインスタンス<1>に対し、時系列データに含まれるデータ数(時刻数)分の残差eを並べた残差ベクトルe<1>と、同個数の信頼度rを並べた信頼度ベクトルr<1>、以下のように得られる。　
　残差ベクトルe<1>=(e1,e2,e3,…..en)　n:はデータの個数
　信頼度ベクトルr<1>=(r1,r2,r3,…,rn)
　なお<>内の数値は、1番目の診断モデルインスタンス<1>に対応することを示す。
　
　次にステップS305cにおいて、信頼度ベクトル内の各要素の値を検査し、所定の信頼閾値以上の値の要素をすべて特定する。そして特定した要素と同じ位置(時刻)に対応する残差を残差ベクトルから抽出する。たとえばrx(xは1～nの任意の値)が信頼閾値以上であれば、残差ベクトルから残差exを抽出する。抽出されなかった残差は信頼性の低いとして以降の処理では用いない。信頼閾値は本発明の第2閾値に対応する。
　次に、センサ情報記憶部21から該当する精度情報を読み出し、読み出した精度情報、およびステップS305cで抽出した残差群を利用して、異常スコア(判定スコア)を求める(S305d)。読み出す精度情報は、たとえば診断モデルインスタンス<1>であれば、風量センサAF1,AF2の精度情報と、VAV1、VAV2の精度情報である。
　図22は異常スコアの計算方法を説明する図である。
　まずステップS305cで抽出した残差群の確率分布Pl（エル）を求める。より詳細に、残差eの平均および分散を計算し、正規分布を仮定して、当該平均および分散に基づいて確率分布Pl（エル）を求める。つまり、Plは、時系列データ記憶部19から読み込んだ時系列データを入力として、l（エル）番目の診断インデックスのパラメータを用いて求められた残差の確率分布である。確率分布Pl（エル）は本発明の第2確率分布に対応する。
　また異常がない場合の確率分布Pn,l（エル）を求める。まず基準時系列データ記憶部に記憶された基準時系列データと、上記最適化したl（エル）番目の診断モデルインスタンス<l>とを用いて、ステップS305b,S305cにより残差群を抽出する。そして、抽出した残差群の平均および分散を求め、正規分布を仮定して、当該平均および分散に基づき確率分布Pn,l（エル）を求める。つまり Pn,lは基準時系列データ17から読み込んだ基準時系列データを入力として、l（エル）番目の診断インデックスのパラメータを用いて求められた残差の確率分布である。確率分布Pn,l（エル）は本発明の第1確率分布に対応する。
　さらに、センサ・VAV(アクチュエータ)毎に静的に異常が発生したと仮定した場合の確率分布Pa,l（エル）,kを求める。kはアクセステーブルの列番号である。本例では、k=1～4であるため、確率分布Pa,l（エル）,1、Pa, l（エル）,2、Pa, l（エル）,3、Pa, l（エル）,4を求める。まず基準時系列データ17の基準時系列データにおける、アクセステーブルのl（エル）行k列の値を、該当する精度情報に応じて補正する。補正した基準時系列データと、上記最適化したl（エル）番目の診断モデルインスタンス<l>とを用いて、ステップS305b,S305cにより残差群を抽出する。そして、抽出した残差分の平均および分散を求め、正規分布を仮定して、当該平均および分散に基づき確率分布Pa,l（エル）,kを求める。これをk=1～4のそれぞれごとに行う。確率分布Pa,l（エル）,kは本発明の第3確率分布に対応する。
　k=1～4のそれぞれに対応するセンサ・アクチュエータの値を補正して残差を求める式を、数式8-1～数式8-4として以下に示す。　
(数式8-1) e(k=1) = (C[i] + α[i, xi +Δxi]) ( Qi / Qj )^2 - (C[j] +α[j, xj])
(数式8-2) e(k=2) = (C[i] + α[i, xi]) ( Qi / Qj )^2 - (C[j] +α[j, xj +Δxj])
(数式8-3) e(k=3) = (C[i] + α[i, xi]) (( Qi + ΔQi )/ Qj )^2 - (C[j] +α[j, xj])
(数式8-4)　e(k=4) = (C[i] + α[i, xi]) ( Qi /(Qj + ΔQj))^2 - (C[j] +α[j, xj])
　基準時系列データにおいて、アクセステーブルのk列の値を補正する様子を図23に示す。(1)が第1列(k=1)の値を補正、(2)が第2列(k=2)、(3)が第3列(k=3)、(4)が第4列(k=4)の値を補正した結果を示す。たとえば(1)では、数式8-1におけるk=1列の指示値xiをxi +Δxiに補正している。補正の方法としては、ここでは、精度情報のプラス方向(あるいはマイナス方向)の最大誤差を加算する。(1)では最大誤差+0．05を加算している。指示値の補正は、開度の最大誤差に相当する指示値の差分を求め、求めた差分を指示値に補正することで行う。(2)～(4)でも同様にして補正を行う。精度情報のプラス方向、マイナス方向の両方を考慮する場合は後述する。
　以上のように各確率分布Pl、Pn,l（エル）、Pa,l（エル）,1、Pa, l（エル）,2、Pa, l（エル）,3、Pa, l（エル）,4を求めたら、拡張KL情報量(Kulback-Leibler divergence)を用いて、異常スコアの計算を行う。
　KL情報量は、確率分布間の差を表す尺度として利用される。KL情報量は対称性がないため距離ではないが、ここで用いる拡張KL情報量は対称性を持つため、確率分布間の距離として定義することができる。
　離散確率変数の確率分布P,Qに関する拡張KL情報量は以下の数式9により定義される。拡張KL情報量は，確率分布P,Q が一致する場合0，ずれが大きいほど値が大きくなり、負の値をとらないことが知られている。
Figure JPOXMLDOC01-appb-M000004
　このKL情報量に基づき、l（エル）番目の最適化された診断モデルインスタンスにおいて、 k 番目のセンサ・アクチュエータに関する異常スコアScore(l,k) を以下のように求める。
　まず、l（エル）番目の最適化された診断モデルインスタンスにおいて、確率分布Pl（エル）と、確率分布Pn,lのずれ(拡張KL情報量)を、下記の数式10により求める。
Figure JPOXMLDOC01-appb-M000005
　またl（エル）番目の最適化された診断インデックスモデルにおいて、確率分布P l（エル）と、確率分布Pa,l（エル）,kのずれ(拡張KL情報量)を、下記の数式11により求める。
Figure JPOXMLDOC01-appb-M000006
　以上に基づき、k 番目のセンサに関する異常スコアScore(l,k) は，以下のように定義することができる。
Figure JPOXMLDOC01-appb-M000007
　図25の上に、診断インデックス<1>～<6>毎に、k=1～4のそれぞれに対応するセンサ・アクチュエータに対して計算された異常スコアを示す。
　上記説明では、補正の際、精度情報のプラス方向およびマイナス方向の一方のみ考慮したが、プラスとマイナス方向の両方を考慮する場合は、以下のようにすればよい。まず、k=1～4に対応する数式8-1～8-4毎に、プラス方向の最大誤差分の補正を行い、計算結果(e)の二乗(e^2)を求める。同様にして、マイナス方向の最大誤差分の補正も行い、計算結果(e)の二乗(e^2)を求める。k=1～4毎に、プラス方向の補正とマイナス方向の補正で得たe^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布Pa,l（エル）,kを求める。一方、基準時系列データに対応する確率分布Pn,l（エル）も、残差eの二乗(e^2)を求め、e^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布を求める。また時系列データ記憶部19内の時系列データに対応する確率分布Pl（エル）も、残差eの二乗(e^2)を求め、e^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布を求める。以降は先と同様にして、拡張KL情報量を求め、異常スコアを計算する。
　次に図18の総合異常判定処理(S306)に進む。
　図24は総合異常判定部134による総合異常判定処理の流れを示すフローチャートである。
　まずセンサ・アクチュエータ(VAV1～VAV3、風量センサAF1～AF3)ごとに、総合スコアを計算する(S306a)。総合スコアは、下記の数式13を用いて計算する。センサ・アクチュエータ毎に、各l（エル）で得られた異常スコアのうち最小値を総合スコアとする。
Figure JPOXMLDOC01-appb-M000008
　図25には総合スコアの計算例が示される。VAV1の開度指示値(ID133)の総合スコアは、診断インデックス<1>～<4>で算出された異常スコア0．96，0.95,0.08,0.09のうちの最小値0．08となる。同様にして他のセンサ・アクチュエータ(ID134,135,143,144,145)についても異常スコアの最小値を選択することにより総合スコアを求める。
　次にステップS306bでは、各センサ・アクチュエータの総合スコアを所定の判定閾値と比較し、総合スコアが判定閾値より大きいときは、該当するセンサ・アクチュエータに対して異常ラベルを設定し(S306c)、判定閾値以下のときは正常ラベルを設定する(S306d)。本例では総合スコアが0.92のVAV2に対して異常ラベルが設定され、他のVAV1,3および風量センサAF1～AF3には正常ラベルが設定される。判定閾値は、本発明の第1閾値に対応する。
　ステップS306eでは、異常検出部13の診断結果として、診断対象期間(時系列データを用いた期間)、計算された異常スコア、総合スコア、センサ・アクチュエータに設定したラベルを、結果格納部20へ書き込む (S306e)。
　上記説明では診断モデルインスタンス<1>～<6>毎に、センサ・アクチュエータの異常スコアを求め、センサ・アクチュエータ毎に最小の異常スコアを選択することで総合スコアを求めた。そして、総合スコアが判定閾値より大きいセンサ・アクチュエータに異常ラベルを設定した。別の実施形態として、まず診断モデルインスタンス<1>～<6>毎に、センサ・アクチュエータの正常スコア(判定スコア)を、下記の数式12aにより求める。そして、センサ・アクチュエータ毎に最大の正常スコアを選択し、総合スコアが判定閾値より小さいセンサ・アクチュエータ（ポイントオブジェクト名）に、異常ラベルを設定し、総合スコアが判定閾値より大きいセンサ・アクチュエータに正常ラベルを設定してもよい。
Figure JPOXMLDOC01-appb-M000009
　異常出部13の異常検出提示情報生成部135は、時系列データ記憶部19、基準時系列データ記憶部17、結果格納部20、精度情報記憶部21に記憶されている情報を利用して、予め定められた形式に従って結果表示画面を動的に生成し、生成した画面を表示部16に表示する。
　図26、図27，図28は異常検出提示情報生成部136が生成した結果表示の画面例である。
　図26は画面内の「リスト」が選択されたときの表示例を示している。総合スコアの大きい順に、各VAVに関するポイントオブジェクト名と、総合スコアが表示される。ただし図示の例では、ポイントオブジェクト名のビル名の表示は省略されている。図左では、異常ラベルが設定されたセンサ・VAV(アクチュエータ)が備え付けられた部屋が網がけ表示されている。
　図27は画面内の「配置図」が選択されたとき表示例を示す。
　図示の画面は、図26の網がけ表示された部屋のうち12階の部屋Cが選択されたときの表示例を示す。図右では、異常ラベルが設定された、VAV2に関するポイントオブジェクト名「Xビル.12.C.I.VAV.2.開度指示値」が網がけ表示されている。
　図28は、画面内の「グラフ」が選択されたときの表示例を示す。
　図28では12階の部屋Cが選択されたときの表示例が示される。図右の上では、横軸が時刻、縦軸が開度指示値の座標系が示される。この座標系に、VAV1の開度指示値のグラフG101, VAV2の開度指示値のグラフG201, VAV3の開度指示値のグラフG301が表示されている。
　また図右の下では、横軸が時刻、縦軸が風量の座標系が示される。この座標系に、VAV1の風量値(風量センサAF1の値)のグラフG102, VAV2の風量値(風量センサAF2の値)のグラフG202, VAV3の風量値(風量センサAF3の値)のグラフG302が表示されている。
　各グラフG101、G201、G301、G102、G202、G302は、診断に用いた日付情報(本例では2009年9月1日0：00～23：59)を利用して、時系列データ記憶部19からは当該期間の時系列データ(本例では1日分)を読み、読み込んだデータをプロットして得られる。
　異常区間Z1はVAV2の開度異常が検出された区間である。詳細を知るため、右下の「解析結果詳細へのリンク」をユーザが入力手段を用いてクリックすると、図29の解析結果画面が生成および表示される。
　この画面も、図26～図28と同様に、異常検出提示情報生成部136が、時系列データ記憶部19、結果格納部20に記憶されている情報を利用して、動的に生成する。
　本例は、Xビル.12.C.I.VAV.2.開度指示値に異常ラベルが設定(総合スコア0.92)されたた場合である。この場合、まずXビル.12Cのモデル情報を読み込み、総合スコア0.92を出力した診断インデックス<2>(図25参照)を参照する。診断インデックス<2>のACCESSのインデックス番号133,134,143,144に該当する時系列データ(診断対象期間分)を、時系列データ記憶部19と基準時系列データ記憶部17からそれぞれ読み込む。
　画面左上のグラフG101、G201、左下のG102,202は、時系列データ記憶部19から読み込んだデータ(観測値)をプロットしたものであり、図28に示したのと同じである。
　画面左上のグラフG211は、VAV2の推定開度指示値である。この推定開度指示値は、数式6を利用して求める。具体的には、数式6のeを0とおき、パラメータを診断モデルインスタンスのPARAMS情報を利用してパラメータを代入する。そして、異常と判断されている x[i] 以外の変数 x[j], Q[i], Q[j] に、時系列データ記憶部19から読み込んだ時系列データを代入してx[i]を求める。x[i] が推定値(期間データ推定値)となる。
　画面右側には回帰プロットが示される。この回帰プロットは以下のようにして作成する。上記同様に数式6の変数 x[j], Q[i], Q[j]に、基準時系列データ記憶部17から読み込んだ基準時系列データを代入することにより、x[i] を、比較データの推定値として得る。そして、当該読み込んだ基準時系列データのI.VAV2.開度指示値に関する値と、当該比較データの推定値の交点に「□」をプロットする。一方、時系列データ記憶部19から読み込んだ時系列データのI.VAV2.開度指示値と、上記で求めた期間データ推定値の交点に「×」をプロットする。
　このように、正常データ(基準時系列データ)が入力されると、観測値と推定値とが略一致するため横軸値=縦軸値(y=x)に漸近するようにデータがプロットされる。一方、異常データが入力されると、観測値から外れる推定値が出力されるため、y=xから逸脱したところにプロットが行われる。
　上記実施形態では、2つのVAVのペア毎に、目的モデルに対し2つの診断モデルインスタンスを生成した。たとえばVAV1とVAV2のペアに対して診断モデルインスタンス＜1＞、＜2＞を生成した。しかしながら、図2のようなシンプルな事例では、ペア毎に、いずれか1つのみを生成しても、異常検出が可能である。したがって、この場合、VAV1とVAV2のペアに対しては、＜1＞および＜2＞のいずれか任意の一方のみを生成し、VAV1とVAV3のペアに対しては、＜3＞および＜4＞のいずれか任意の一方のみを生成し、VAV2とVAV3のペアに対しては、＜5＞および＜6＞のいずれか任意の一方のみを生成してもよい。
　また上記実施形態では、目的モデルF01に従って、2つのVAVのすべてのペアを作成し、ペア毎に診断モデルインスタンスを生成したが、本発明はこれに限定されない。たとえば多数のVAVが存在するときに、相互の影響が弱いと事前に分かっているペアについては、診断モデルインスタンスを生成しなくてもよい。診断モデルインスタンスを生成すべきペアは事前に指定しておけばよい。目的モデルF01ではすべてのペアを生成することが暗に示されている。
　以上のように、本実施形態によれば、個体差のある未知の特性を持つ機器（センサ・アクチュエータ）の異常検出を高い精度で行うことが可能となる。
　
(変形例)
　これまで説明した実施形態では、ダクトに設置されたVAVの開度指示値と、VAV風量を事例として取り上げたが、本発明は、これに限定されず、たとえば以下のような事例にも適用可能である。
　たとえば水量調整装置における冷水弁の開度指示値(vw)と、冷水弁の水量(Qw)の事例も可能である。この場合、数式(1)から、以下のようなモデルインスタンスを生成できる。　
　　min_{Ci,Cj,βi,βj } 
 R　=　{ (Ci + βi(vwi) ) × (Qwi/Qwj)^2 － (Cj +βj(vwj)) }^2
　以上は一例であり、本発明は他の種々の事例にも適用可能である。

Claims (8)

  1. 　媒体を流通させる媒体経路における異なる箇所に配置され、それぞれ外部から与えられた指示値に従って前記媒体経路内の前記媒体を共通に出力制御する第1～第K制御機器と、前記第1～第K制御機器によって出力制御される媒体の出力量を計測する第1～第Kセンサに対する異常検出装置であって、
    　(Ci + βi(Xi) ) × (Qi/Qj)^2 － (Cj + βj(Xj) (Xiは第1観測変数、Bi(Xi)は Xiを入力変数とする未知の関数である第1未知特性項、Ciは第1定数項、Qiは第2観測変数、Xjは第3観測変数、Bj(Xj)は Xjを入力変数とする未知の関数である第2未知特性項、Cjは第2定数項、Qjは第4観測変数)によって定義される目的モデルを記憶する目的モデル格納部と、
    　第1の期間において前記第1～第K制御機器に与えられた前記指示値と、前記第1～第Kセンサによる計測値とを含む第1時系列データを所定時間間隔で記憶する第1記憶部と、
    　前記第1～第K制御機器における2つの制御機器のペア毎に、前記目的モデルの前記Xiに一方の制御機器の指示値、前記Ciに前記一方の制御機器の配置箇所に応じた未知の定数、前記Qiに前記一方の制御機器に対応するセンサの計測値、前記Xjに他方の制御機器の指示値、前記Cjに前記他方の制御機器の配置箇所に応じた未知の定数、前記Qjに前記他方の制御機器に対応するセンサの計測値を割り当てることによって、それぞれ診断モデルインスタンスを生成し、
    　前記ペア毎に生成したすべての診断モデルインスタンスのそれぞれにおける前記Bi(Xi)およびBj(Xj)を、前記XiおよびXjの値の範囲をそれぞれ分割したセクション毎に一定の値を取る関数と定義し、前記第1記憶部内の前記第1時系列データを用いて、前記すべての診断モデルインスタンスの値の二乗を最小化するように前記Bi(Xi)、前記Ci、前記Bj(Xj)、前記Cjを同定することにより、最適化された診断モデルインスタンスをそれぞれ得るモデル最適化部と、
    　前記第1の期間と異なる第2の期間において、前記第1～第K制御機器に与えられた前記指示値と、前記第1～第Kセンサによる計測値とを含む第2時系列データを所定時間間隔で記憶する第2記憶部と、
    　前記第1～第Kセンサによる最大観測誤差と、前記指示値に対する前記第1～第K制御機器の最大応答誤差とを表す精度情報を記憶する第3記憶部と、
    　前記最適化された診断モデルインスタンス毎に、
    　前記第1記憶部内の前記第1時系列データを用いて前記最適化された診断モデルインスタンスの値である第1残差を求め、前記第1残差の確率分布である第1確率分布を生成し、
    　前記第2記憶部内の前記第2時系列データを用いて前記最適化された診断モデルインスタンスの値である第2残差を求め、前記第2残差の確率分布である第2確率分布を生成し、
    　前記最適化された診断モデルインスタンスを生成する元となった2つの制御機器および2つのセンサのそれぞれ毎に、前記第1時系列データにおけるそれぞれに対応する値をそれぞれの前記精度情報に応じて変更し、変更後の時系列データを用いて前記最適化された診断モデルインスタンスの値である第3残差を求め、前記第3残差の確率分布である第3確率分布をそれぞれ生成し、
    　前記第2確率分布と前記第1確率分布間の距離を表す第1の拡張KL(Kulback-Leibler divergence)情報量と、前記第2確率分布と各前記第3確率分布のそれぞれとの距離を表す第2の拡張KL情報量とを求め、前記第1および第2拡張KL情報量の合計に対する前記第1拡張KL情報量の比率に応じて判定スコアをそれぞれ計算する、
    　判定スコア計算部と、
    　前記第1～第K制御機器および前記第1～Kセンサ毎に、前記最適化された診断モデルインスタンスから得られた判定スコアの最小値である総合スコアを求め、第1閾値より大きい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する総合異常判定部と、
    　を備えた異常検出装置。
  2. 　前記モデル最適化部は、前記第1～第K制御機器毎に、前記セクションに前記指示値が含まれるデータ数をカウントし、データ数が多いほど大きな信頼値を前記セクションに設定し、
    　前記判定スコア計算部は、
    　前記第1時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、第2閾値以上の前記信頼度が計算されたデータから求められた第1残差のみを用いて前記第1確率分布を生成し、
    　前記第2時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第2残差のみを用いて前記第2確率分布を生成し、
    　前記変更後の時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第3残差のみを用いて前記第3確率分布を生成する、
    　ことを特徴とする請求項1に記載の異常検出装置。
  3. 　前記モデル最適化部は、さらに、前記Xiに他方の制御機器の指示値、前記Ciに前記他方の制御機器に関する未知の定数、前記Qiに前記他方の制御機器に対応するセンサの計測値、前記Xjに前記一方の制御機器の指示値、前記Cjに前記一方の制御機器に関する未知の定数、前記Qjに前記一方の制御機器に対応するセンサの計測値を割り当てることによって診断モデルインスタンスを生成する
    　ことを特徴とする請求項1に記載の異常検出装置。
  4. 　前記判定スコア計算部は、
    　前記第1残差の二乗を求め、前記第1残差の二乗の確率分布を前記第1確率分布として生成し、
    　前記第2残差の二乗を求め、前記第2残差の二乗の確率分布を第2確率分布として生成し、
    　前記第3残差の二乗を求め、前記第3残差の二乗の確率分布である第3確率分布を生成する
    　ことを特徴とする請求項1に記載の異常検出装置。
  5. 　前記判定スコア計算部は、前記第1および第2拡張KL情報量の合計に対する前記第2拡張KL情報量の比率に応じて前記判定スコアを計算し、
    　前記総合異常判定部は、前記判定スコアの最大値を前記総合スコアとして求め、前記第1閾値より小さい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する
    　ことを特徴とする請求項1に記載の異常検出装置。
  6. 　前記目的モデル格納部は、前記目的モデルとして、前記(Ci + βi(Xi) ) × (Qi/Qj)^2 － (Cj + βj(Xj)に代えて、(Ci + βi(Xi)) × f(Qi,Qj) － (Cj + βj(Xj) ) (f(Qi,Qj)は、前記Qi,Qjを入力変数とする所定の関数)を記憶する
    　ことを特徴とする請求項1に記載の異常検出装置。
  7. 　前記第1～第K制御機器は、可変風量装置または固定風量装置であり、前記第1～第K制御機器に与えられる指示値は、前記可変風量装置または前記固定風量装置の弁の開度指示値であり、前記第1～Kセンサは風量センサであることを特徴とする請求項1に記載の異常検出装置。
  8. 　前記第1～第K制御機器は、水量調整装置であり、前記第1～第K制御機器に与えられる指示値は、前記水量調整装置の弁の開度指示値であり、前記第1～Kセンサは水量センサであることを特徴とする請求項1に記載の異常検出装置。
PCT/JP2010/055695 2010-03-30 2010-03-30 異常検出装置 WO2011121726A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/JP2010/055695 WO2011121726A1 (ja) 2010-03-30 2010-03-30 異常検出装置
JP2012507951A JP5337909B2 (ja) 2010-03-30 2010-03-30 異常検出装置
US13/628,951 US8577649B2 (en) 2010-03-30 2012-09-27 Anomaly detecting apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2010/055695 WO2011121726A1 (ja) 2010-03-30 2010-03-30 異常検出装置

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US13/628,951 Continuation US8577649B2 (en) 2010-03-30 2012-09-27 Anomaly detecting apparatus

Publications (1)

Publication Number Publication Date
WO2011121726A1 true WO2011121726A1 (ja) 2011-10-06

Family

ID=44711517

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2010/055695 WO2011121726A1 (ja) 2010-03-30 2010-03-30 異常検出装置

Country Status (3)

Country Link
US (1) US8577649B2 (ja)
JP (1) JP5337909B2 (ja)
WO (1) WO2011121726A1 (ja)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015079975A1 (ja) * 2013-11-26 2015-06-04 株式会社トクヤマ 運動装置の状態監視システム
JP2015179440A (ja) * 2014-03-19 2015-10-08 株式会社東芝 センサ割当て装置及びセンサ診断装置
JP2017021790A (ja) * 2015-06-22 2017-01-26 ジーイー・アビエイション・システムズ・リミテッドＧｅ Ａｖｉａｔｉｏｎ Ｓｙｓｔｅｍｓ Ｌｉｍｉｔｅｄ 隠れマルコフモデルの混合を使用する検証および異常検出のためのシステムおよび方法
JP2017083985A (ja) * 2015-10-26 2017-05-18 株式会社Ｓｃｒｅｅｎホールディングス 時系列データ処理方法、時系列データ処理プログラム、および、時系列データ処理装置
CN109997087A (zh) * 2016-12-01 2019-07-09 住友重机械工业株式会社 故障诊断系统
JP2020035458A (ja) * 2016-02-05 2020-03-05 株式会社東芝 情報処理装置、方法、及びプログラム
WO2022065309A1 (ja) * 2020-09-23 2022-03-31 ダイキン工業株式会社 情報処理装置、情報処理方法、及びプログラム

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9197511B2 (en) * 2012-10-12 2015-11-24 Adobe Systems Incorporated Anomaly detection in network-site metrics using predictive modeling
FR3010448B1 (fr) 2013-09-06 2015-08-21 Snecma Procede de surveillance d’une degradation d’un dispositif embarque d’un aeronef avec determination automatique d’un seuil de decision
IL229819A (en) * 2013-12-05 2016-04-21 Deutsche Telekom Ag System and method for anomalously identifying information technology servers through incident consolidation
US20150247584A1 (en) * 2014-03-03 2015-09-03 Brent Richard SINGLEY Flood prevention device
JP2015179443A (ja) * 2014-03-19 2015-10-08 株式会社東芝 診断モデル生成装置、診断用モデル生成方法、及び異常診断装置
US9727533B2 (en) * 2014-05-20 2017-08-08 Facebook, Inc. Detecting anomalies in a time series
US9680646B2 (en) * 2015-02-05 2017-06-13 Apple Inc. Relay service for communication between controllers and accessories
JP6350736B2 (ja) 2015-02-17 2018-07-04 富士通株式会社 判定装置、判定方法および判定プログラム
US9823160B2 (en) * 2015-04-02 2017-11-21 The Boeing Company Apparatus and methods for testing suction cups mounted to a track
KR20170007958A (ko) * 2015-07-13 2017-01-23 에스케이하이닉스 주식회사 메모리 시스템
JP6603078B2 (ja) 2015-08-31 2019-11-06 Ｎｔｎ株式会社 車輪用軸受装置
US10685306B2 (en) * 2015-12-07 2020-06-16 Sap Se Advisor generating multi-representations of time series data
US10095757B2 (en) 2015-12-07 2018-10-09 Sap Se Multi-representation storage of time series data
US9785495B1 (en) * 2015-12-14 2017-10-10 Amazon Technologies, Inc. Techniques and systems for detecting anomalous operational data
US10642813B1 (en) 2015-12-14 2020-05-05 Amazon Technologies, Inc. Techniques and systems for storage and processing of operational data
US11037066B2 (en) 2016-07-13 2021-06-15 International Business Machines Corporation Estimation of abnormal sensors
US10628252B2 (en) * 2017-11-17 2020-04-21 Google Llc Real-time anomaly detection and correlation of time-series data
JP6928669B2 (ja) * 2017-11-28 2021-09-01 株式会社安川電機 制御システム、工場システム、学習システム、推定用モデルの生成方法及びアクチュエータの状態推定方法
US11036715B2 (en) * 2018-01-29 2021-06-15 Microsoft Technology Licensing, Llc Combination of techniques to detect anomalies in multi-dimensional time series
JP7038629B2 (ja) * 2018-08-31 2022-03-18 三菱電機ビルテクノサービス株式会社 機器状態監視装置及びプログラム
US11792215B1 (en) * 2018-12-04 2023-10-17 Amazon Technologies, Inc. Anomaly detection system for a data monitoring service
EP3796117B1 (de) * 2019-09-18 2021-10-27 Siemens Aktiengesellschaft Diagnoseverfahren und diagnosesystem für eine verfahrenstechnische anlage
CN111459778B (zh) * 2020-03-12 2024-05-07 平安科技（深圳）有限公司 运维系统异常指标检测模型优化方法、装置及存储介质
CN112905671A (zh) * 2021-03-24 2021-06-04 北京必示科技有限公司 时间序列异常处理方法、装置、电子设备及存储介质
US20230259113A1 (en) * 2022-02-11 2023-08-17 Novity, Inc. Subsystem-level model-based diagnosis
US11821324B2 (en) 2022-04-25 2023-11-21 General Electric Company Duct failure detection in a turbine engine

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003131733A (ja) * 2001-10-26 2003-05-09 Kajima Corp 設備制御監視方法及び設備制御監視装置
JP2003208219A (ja) * 2002-01-10 2003-07-25 Yamatake Building Systems Co Ltd 異常検知装置および方法
JP2008065393A (ja) * 2006-09-04 2008-03-21 Research Organization Of Information & Systems グループ判別装置及びグループ判別方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3254624B2 (ja) 1996-05-31 2002-02-12 株式会社山武 スティックスリップ検出方法および検出装置
US6223544B1 (en) 1999-08-05 2001-05-01 Johnson Controls Technology Co. Integrated control and fault detection of HVAC equipment
US7590513B2 (en) * 2006-01-30 2009-09-15 Nec Laboratories America, Inc. Automated modeling and tracking of transaction flow dynamics for fault detection in complex systems
JP4413915B2 (ja) * 2006-12-13 2010-02-10 株式会社東芝 異常兆候検出装置および方法
JP4782727B2 (ja) * 2007-05-17 2011-09-28 株式会社東芝 機器状態監視装置並びに機器状態監視のための方法およびプログラム
WO2010082322A1 (ja) * 2009-01-14 2010-07-22 株式会社日立製作所 装置異常監視方法及びシステム

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003131733A (ja) * 2001-10-26 2003-05-09 Kajima Corp 設備制御監視方法及び設備制御監視装置
JP2003208219A (ja) * 2002-01-10 2003-07-25 Yamatake Building Systems Co Ltd 異常検知装置および方法
JP2008065393A (ja) * 2006-09-04 2008-03-21 Research Organization Of Information & Systems グループ判別装置及びグループ判別方法

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10428813B2 (en) 2013-11-26 2019-10-01 Tokuyama Corporation State monitoring system of moving apparatus
JPWO2015079975A1 (ja) * 2013-11-26 2017-03-16 株式会社トクヤマ 運動装置の状態監視システム
WO2015079975A1 (ja) * 2013-11-26 2015-06-04 株式会社トクヤマ 運動装置の状態監視システム
US10627265B2 (en) 2014-03-19 2020-04-21 Kabushiki Kaisha Toshiba Sensor assignment apparatus and sensor diagnostic apparatus
JP2015179440A (ja) * 2014-03-19 2015-10-08 株式会社東芝 センサ割当て装置及びセンサ診断装置
JP2017021790A (ja) * 2015-06-22 2017-01-26 ジーイー・アビエイション・システムズ・リミテッドＧｅ Ａｖｉａｔｉｏｎ Ｓｙｓｔｅｍｓ Ｌｉｍｉｔｅｄ 隠れマルコフモデルの混合を使用する検証および異常検出のためのシステムおよび方法
JP2017083985A (ja) * 2015-10-26 2017-05-18 株式会社Ｓｃｒｅｅｎホールディングス 時系列データ処理方法、時系列データ処理プログラム、および、時系列データ処理装置
JP2020035458A (ja) * 2016-02-05 2020-03-05 株式会社東芝 情報処理装置、方法、及びプログラム
CN109997087A (zh) * 2016-12-01 2019-07-09 住友重机械工业株式会社 故障诊断系统
US11422544B2 (en) 2016-12-01 2022-08-23 Sumitomo Heavy Industries, Ltd. Failure diagnosis system
JP2022052545A (ja) * 2020-09-23 2022-04-04 ダイキン工業株式会社 情報処理装置、情報処理方法、及びプログラム
WO2022065309A1 (ja) * 2020-09-23 2022-03-31 ダイキン工業株式会社 情報処理装置、情報処理方法、及びプログラム
JP7093031B2 (ja) 2020-09-23 2022-06-29 ダイキン工業株式会社 情報処理装置、情報処理方法、及びプログラム
US11774955B1 (en) 2020-09-23 2023-10-03 Daikin Industries, Ltd. Information processing apparatus, information processing method, and program

Also Published As

Publication number Publication date
US8577649B2 (en) 2013-11-05
JPWO2011121726A1 (ja) 2013-07-04
JP5337909B2 (ja) 2013-11-06
US20130024172A1 (en) 2013-01-24

Similar Documents

Publication Publication Date Title
JP5337909B2 (ja) 異常検出装置
JP6095732B2 (ja) コンピュータベースの方法及びコンピュータベースのモデル
US7853431B2 (en) On-line monitoring and diagnostics of a process using multivariate statistical analysis
EP2791745B1 (en) A method of operating a process or machine
WO2008085705A1 (en) Method and system for modeling a process in a process plant
JP6240050B2 (ja) バイアス推定装置、その方法およびプログラム、並びに故障診断装置、その方法およびプログラム
WO2008014344A2 (en) Model based method for detecting abnormal operation of a level regulatory control loop and associated apparatus
Wang et al. Fault detection and diagnosis for multiple faults of VAV terminals using self-adaptive model and layered random forest
CN113287104A (zh) 数据分类装置
JP5264796B2 (ja) プラント運転支援装置
AU2024201310A1 (en) Methods for identifying key performance indicators
Qi et al. A Bayesian approach for control loop diagnosis with missing data
JP6448745B2 (ja) バイアス推定装置、その方法およびプログラム
Rusinov et al. Fault diagnosis in chemical processes with application of hierarchical neural networks
JP2010146137A (ja) パラメータ調整支援装置
US11669082B2 (en) Online fault localization in industrial processes without utilizing a dynamic system model
Blázquez et al. Neuro-fuzzy identification applied to fault detection in nonlinear systems
Lin et al. Implementation and test of an automated control hunting fault correction algorithm in a fault detection and diagnostics tool
Pradhan A Dynamic Bayesian Network Framework for Data-Driven Fault Diagnosis and Prognosis of Smart Building Systems
Dedemen et al. Quantifying performance degradation of HVAC systems for proactive maintenance using a data-driven approach
WO2008042757A2 (en) Univariate method for monitoring and analysis of multivariate data
JP2022103931A (ja) モデルを生成する方法、プログラムおよび装置
JP2023176898A (ja) ニューラルネットワークの計算装置及び計算方法
CN117043700A (zh) 基于模拟的征兆的异常的原因分析
Rato et al. A systematic comparison of statistical process monitoring methods for high-dimensional, time-dependent processes

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 10848905

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2012507951

Country of ref document: JP

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 10848905

Country of ref document: EP

Kind code of ref document: A1