WO2022130611A1 - 情報処理装置及び情報処理方法 - Google Patents

情報処理装置及び情報処理方法 Download PDF

Info

Publication number
WO2022130611A1
WO2022130611A1 PCT/JP2020/047394 JP2020047394W WO2022130611A1 WO 2022130611 A1 WO2022130611 A1 WO 2022130611A1 JP 2020047394 W JP2020047394 W JP 2020047394W WO 2022130611 A1 WO2022130611 A1 WO 2022130611A1
Authority
WO
WIPO (PCT)
Prior art keywords
time point
sensor data
value
time
unit
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/JP2020/047394
Other languages
English (en)
French (fr)
Inventor
直人 高野
洋一 堀澤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Mitsubishi Electric Corp
Original Assignee
Mitsubishi Electric Corp
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 Mitsubishi Electric Corp filed Critical Mitsubishi Electric Corp
Priority to US18/031,172 priority Critical patent/US12498709B2/en
Priority to KR1020237018898A priority patent/KR102841817B1/ko
Priority to JP2021526355A priority patent/JP6935046B1/ja
Priority to DE112020007851.5T priority patent/DE112020007851T5/de
Priority to PCT/JP2020/047394 priority patent/WO2022130611A1/ja
Priority to CN202080107328.7A priority patent/CN116569120A/zh
Priority to TW110146274A priority patent/TWI814170B/zh
Publication of WO2022130611A1 publication Critical patent/WO2022130611A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M15/00—Testing of engines
    • G01M15/04—Testing internal-combustion engines
    • G01M15/042—Testing internal-combustion engines by monitoring a single specific parameter not covered by groups G01M15/06 - G01M15/12
    • G01M15/046—Testing internal-combustion engines by monitoring a single specific parameter not covered by groups G01M15/06 - G01M15/12 by monitoring revolutions
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00—Subject matter not provided for in other groups of this subclass
    • 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/024—Quantitative history assessment, e.g. mathematical relationships between available data; Functions therefor; Principal component analysis [PCA]; Partial least square [PLS]; Statistical classifiers, e.g. Bayesian networks, linear regression or correlation analysis; 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
    • 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/0245—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 qualitative model, e.g. rule based; if-then decisions

Definitions

  • This disclosure relates to an information processing device and an information processing method for performing information processing for diagnosing the state of a mechanical device.
  • a mechanical device including a mechanical part such as a ball screw or a speed reducer
  • the mechanical part deteriorates over time, and various abnormalities such as increased friction, vibration, and damage to the housing occur. Therefore, an information processing device that detects and grasps this kind of abnormality at an early stage is considered to be important for efficient factory operation.
  • Patent Document 1 An example of an information processing device for diagnosing the state of a mechanical device is described in Patent Document 1 below.
  • Patent Document 1 describes a method of holding sensor data acquired in time series for a preset period and calculating feature quantities such as mean value and variance from the held sensor data. Further, in Patent Document 1, the calculated feature amount is retained for a preset period, an index such as skewness is calculated based on the retained feature amount, and the calculated value is compared with the past value. It is described that the abnormality of the mechanical device is detected in.
  • the present disclosure has been made in view of the above, and obtains an information processing apparatus capable of performing information processing for diagnosing the state of a mechanical device without using a computer having a high computing power and a computer having a large storage capacity.
  • the purpose is.
  • the information processing apparatus includes a sensor data acquisition unit, an internal variable holding unit, an internal variable calculation unit, a feature quantity calculation unit, and a state diagnosis unit.
  • the sensor data acquisition unit acquires the measured value of the physical quantity of the mechanical device measured by the sensor, and the sensor data is the measured value from the first time point to the Nth time point (N is an integer of 2 or more) among the measured values.
  • the internal variable holding unit holds a number of internal variables smaller than N, which are sequentially calculated in time series based on the sensor data.
  • the internal variable calculation unit calculates the internal variable corresponding to the j + 1th time point (j is an integer from 1 to N-1) based on the sensor data at the j + 1st time point and the internal variable corresponding to the jth time point.
  • the feature amount calculation unit calculates the feature amount obtained by extracting the statistical features included in the sensor data from the first time point to the Nth time point based on the internal variables of the Nth time point.
  • the state diagnosis unit diagnoses the state of the mechanical device based on the feature amount.
  • FIG. 1 A flowchart used to explain the information processing method according to the first embodiment.
  • FIG. 1 is a diagram showing a configuration example of an information processing system 100 including an information processing device 1000 according to the first embodiment.
  • the information processing system 100 includes an information processing device 1000, a mechanical device 1008, a sensor 1010, a motor 1009, and a device control unit 1099.
  • the motor 1009 drives the mechanical device 1008 by applying a driving force to the mechanical device 1008.
  • the sensor 1010 measures the physical quantity of the mechanical device 1008. Examples of physical quantities measured by the mechanical device 1008 are position, velocity, acceleration, operation command, current, voltage, torque, force, pressure, voice, or light quantity.
  • the motor speed and the motor torque will be described as an example.
  • the motor speed is the rotation speed of the motor 1009
  • the motor torque is the torque generated by the motor 1009.
  • the sensor 1010 outputs a sensor signal including a measured value of a physical quantity to the information processing device 1000 and the device control unit 1099.
  • the device control unit 1099 determines a control signal for controlling the motor 1009 based on the sensor signal.
  • the motor 1009 is controlled by a control signal output from the device control unit 1099.
  • the information processing device 1000 includes a sensor data acquisition unit 1001, an internal variable holding unit 1002, an internal variable calculation unit 1003, a feature quantity calculation unit 1004, an initialization processing unit 1005, and a state diagnosis unit 1006.
  • the sensor data acquisition unit 1001 acquires the measured value of the physical quantity of the mechanical device 1008 measured by the sensor 1010. As described above, the measured value of the physical quantity is included in the sensor signal transmitted by the sensor 1010. Further, the sensor data acquisition unit 1001 holds the measured values from the first time point to the Nth time point among the acquired measured values as sensor data.
  • N is an integer of 2 or more.
  • the internal variable holding unit 1002 temporarily holds a smaller number of internal variables than N. Internal variables are variables that are sequentially calculated in time series based on sensor data. Internal variables are used to calculate features. Details of internal variables and features will be described later.
  • the internal variable calculation unit 1003 receives the sensor data transmitted from the sensor data acquisition unit 1001 and the internal variables transmitted from the internal variable holding unit 1002.
  • the internal variable calculation unit 1003 updates the internal variable based on the sensor data and the internal variable. More generalized, the internal variable calculation unit 1003 calculates the internal variable corresponding to the j + 1th time point based on the sensor data at the j + 1th time point and the internal variable corresponding to the jth time point, so that the internal variable can be calculated. Will be updated sequentially.
  • j is an integer from 1 to N-1. That is, the internal variable calculation unit 1003 calculates the internal variable at a certain time point based on the sensor data at a certain time point and the internal variable one time point before the certain time point.
  • the calculated internal variable is transmitted to the internal variable holding unit 1002.
  • the feature amount calculation unit 1004 receives the sensor data transmitted from the sensor data acquisition unit 1001 and the internal variables transmitted from the internal variable holding unit 1002. The feature amount calculation unit 1004 calculates the feature amount based on the sensor data and the internal variables. More generally, the feature amount calculation unit 1004 calculates the feature amount obtained by extracting the statistical features included in the sensor data from the first time point to the Nth time point based on the internal variables. The calculated feature amount is transmitted to the state diagnosis unit 1006.
  • the initialization processing unit 1005 executes the initialization processing.
  • the initialization process is a process of setting an internal variable held by the internal variable holding unit 1002 to an initial value. More generally, the initialization processing unit 1005 performs a process of determining the internal variable at the first time point to a value between the preset maximum value and the preset minimum value.
  • the state diagnosis unit 1006 performs a diagnostic process for diagnosing the state of the mechanical device 1008 based on the feature amount, and outputs a diagnostic result which is the result of the diagnostic process.
  • FIG. 2 is a diagram showing a hardware configuration example of the mechanical device 1008 and its peripheral device according to the first embodiment.
  • FIG. 2 shows a mechanical device 1008 using a servomotor 1230 as a drive source as a configuration example of the mechanical device 1008.
  • the drive torque generated by the servo motor 1230 is output from the servo motor shaft 1231 and input to the ball screw shaft 1224 via the coupling 1220.
  • the ball screw 1210 converts the rotational movement into a linear movement by a screw mechanism, and operates the movable portion 1212.
  • the movable portion 1212 is connected to different mechanical parts, and the moved mechanical parts are used according to the purpose of the mechanical device 1008.
  • the movable portion 1212 is restricted from moving in a desired direction by the guide 1213.
  • the guide 1213 assists the movable portion 1212 so that the mechanical device 1008 can operate with high accuracy.
  • the servomotor 1230 includes an encoder 1233 that measures the rotation angle and a current sensor 1232 that measures the current so that the servomotor shaft 1231 can be driven following a predetermined position, speed, or torque. It is common that and is attached.
  • the driver 1240 performs feedback control based on the information obtained from the current sensor 1232, and supplies the electric power required for driving to the servomotor 1230. The calculation required for feedback control is performed by the device control unit 1099 in FIG.
  • the motor torque is exemplified as an example of the sensor data used for diagnosing the state, but the present invention is not limited to this. Any physical quantity such as position, velocity, acceleration, current, voltage, torque, force, pressure, voice, and light amount described above may be used as long as the state of the mechanical device 1008 is included as information. Further, instead of these physical quantities, image information or the like may be used.
  • the current sensor 1232 and the encoder 1233 are illustrated, but the present invention is not limited thereto.
  • a laser displacement sensor, a gyro sensor, a vibrometer, an acceleration sensor, a voltmeter, a torque sensor, a pressure sensor, a microphone, an optical sensor, a camera and the like can be exemplified.
  • the mounting position of the sensor 1010 does not necessarily have to be close to the servomotor 1230, and may be a position suitable for diagnosing the state of the mechanical device 1008.
  • an acceleration sensor may be installed on the outer surface of the guide 1213 or the like, and the acceleration may be measured as sensor data.
  • the PLC (Programmable Logic Controller) 1260 sends an operation command of the servomotor 1230 to the driver 1240.
  • a PC Personal Computer
  • PC1270 may be prepared as needed.
  • PC1270 is used to send a command to PLC1260.
  • an industrial PC Fractory Automation PC or Industrial PC
  • a PLC display 1250 for monitoring the status of the PLC 1260 and a PC display 1280 for monitoring the status of the PC 1270 may be prepared.
  • a plurality of drive sources such as the servomotor 1230 are provided in one mechanical device 1008. Therefore, a plurality of drivers 1240 may be prepared as needed.
  • a single PLC1260 may control or a plurality of PLC1260s may cooperate to operate the mechanical device 1008. Also in the case of these configurations, the information processing apparatus 1000 shown in the present embodiment can be implemented in the same manner.
  • FIG. 3 is a diagram showing a configuration example in the case where the processing circuit included in the driver 1240 shown in FIG. 2 is configured by the processor 1291 and the memory 1292.
  • the processing circuit is composed of a processor 1291 and a memory 1292
  • each function of the processing circuit of the driver 1240 is realized by software, firmware, or a combination of software and firmware.
  • the software or firmware is written as a program and stored in memory 1292.
  • each function is realized by the processor 1291 reading and executing the program stored in the memory 1292. That is, the processing circuit includes a memory 1292 for storing a program in which the processing of the driver 1240 will be executed as a result. It can also be said that these programs cause the computer to execute the procedures and methods of the driver 1240.
  • the processor 1291 may be a computing means called a CPU (Central Processing Unit), a processing device, a computing device, a microprocessor, a microcomputer, or a DSP (Digital Signal Processor).
  • the memory 1292 may be, for example, a non-volatile or volatile semiconductor memory such as a RAM, a ROM (Read Only Memory), a flash memory, an EPROM (Erasable Project ROM), or an EEPROM (registered trademark) (Electrically EPROM). Further, the memory 1292 may be used as a storage means such as a magnetic disk, a flexible disk, an optical disk, a compact disk, a mini disk, or a DVD (Digital Versaille Disc).
  • FIG. 4 is a diagram showing a configuration example in the case where the processing circuit included in the driver 1240 shown in FIG. 2 is configured with dedicated hardware.
  • the processing circuit 1293 shown in FIG. 4 may be, for example, a single circuit, a composite circuit, a programmed processor, a parallel programmed processor, an ASIC (Application Specific Integrated Circuit), or the like. FPGA (Field Processor Metal Gate Array) or a combination thereof may be used.
  • the functions of the driver 1240 may be realized by the processing circuit 1293 for each function, or a plurality of functions may be collectively realized by the processing circuit 1293.
  • the driver 1240 and the PLC1260 may be connected via a network. Further, the PC1270 may exist on the cloud server.
  • An example of the hardware configuration is as described above, but the driver 1240, PLC1260, and PC1270 are not indispensable, and a device for implementing the information processing device according to the present disclosure is separately prepared and implemented inside the device. May be good.
  • a single device having a battery, a microcomputer, a sensor, a display, and a communication function may be used, and the state of the mechanical device 1008 may be estimated based on the sensor data obtained by acquiring the sound of the mechanical device 1008 with a microphone. good.
  • the indicators are not essential, and instead of displaying the results on the PLC display 1250 or the PC display 1280, the driver 1240, the existing LED provided in the PLC 1260, etc. are used to display the results. It may be displayed. Further, the servomotor 1230 may be stopped driving when it is diagnosed that an abnormality has occurred without showing the result on the display.
  • the motor speed is obtained from the encoder 1233 provided in the servomotor 1230, but the present invention is not limited to this.
  • the motor speed may be obtained by using a control signal from the device control unit 1099 that gives a drive command to the motor 1009.
  • the servo motor 1230 has been described as a rotary servo motor, but other motors or drive sources such as linear servo motors, inducer motors, stepping motors, brush motors, ultrasonic motors, etc. May be carried out using.
  • the ball screw 1210 and the coupling 1220 are examples of components and are not limited thereto.
  • the information processing device according to the present disclosure can also be applied to a mechanical device composed of various other parts such as a speed reducer, a guide, a belt, a screw, a pump, a bearing, and a housing.
  • FIG. 5 is a diagram showing time-series waveforms of the motor speed and the motor torque in the first embodiment.
  • the signal directly obtained from the current sensor 1232 is a signal obtained by measuring the three-phase current flowing in the motor 1009.
  • the motor torque can be calculated by applying appropriate conversion to the three-phase current.
  • the three-phase current values detected by the current sensor 1232 may be used for diagnosis as sensor data.
  • the signal obtained from the encoder 1233 is position information representing the rotation angle of the motor 1009. Therefore, if the position information is subjected to processing such as numerical differentiation, the motor speed, which is the rotation speed of the motor 1009, can be obtained. Therefore, instead of the signal obtained from the current sensor 1232, the signal obtained from the encoder 1233 may be used for diagnosis as sensor data.
  • FIGS. 5 (a) and 5 (b) represents time.
  • FIG. 5A shows a time-series waveform of the motor speed
  • FIG. 5B shows a time-series waveform of the motor torque. These are the waveforms when a single drive called positioning is performed from the state where the motor 1009 is stopped.
  • the sampling cycle which is the acquisition cycle of the waveform in FIG. 5B, is 0.5 milliseconds
  • the sampling period, acquisition time, and data score N shown here are examples, and are not limited to these numerical values.
  • the number of positionings is one, but the number of positionings may be a plurality of times. Further, although the operation when positioning the motor 1009 will be described here, the present invention is not limited to this. This embodiment can also be applied to controls other than positioning, such as speed control or torque control.
  • the motor speed of FIG. 5A will be described.
  • the motor 1009 is stopped between the time Tr0 and the time Tr1, and the motor speed is 0 [r / min]. This period is referred to as "Ts1”.
  • the motor 1009 accelerates between the time Tr1 and the time Tr2, and the motor speed increases to 500 [r / min]. This period is referred to as "Ta”.
  • the motor speed is constant and remains at 500 [r / min]. This period is referred to as "Te”.
  • the motor 1009 decelerates between the time Tr3 and the time Tr4, and the motor speed decreases to 0 [r / min].
  • the motor torque shown in FIG. 5 (b) will be described.
  • the motor 1009 is stopped, and the motor torque required for the operation of the motor 1009 is almost 0 [Nm].
  • a torque for accelerating the mechanical device 1008 is required, and a motor torque larger than the period Ts1 is generated. If viscous friction is present in the mechanical device 1008 due to the influence of viscous friction or the like, the required torque may gradually increase depending on the speed, as shown in FIG. In the following period Te, the speed of the mechanical device 1008 does not change, and a substantially constant motor torque is generated.
  • the motor torque described above is the torque required for the operation of the almost ideal mechanical device 1008.
  • a large or small noise as shown in FIG. 5B is generated in the motor torque.
  • the ball screw shaft 1224, guide 1213, coupling 1220, etc. which are parts of the mechanical device 1008, deteriorate, the magnitude of the above-mentioned friction or vibration changes, and the influence of friction or vibration appears on the sensor data including the motor torque. It is known. Therefore, the sensor data is analyzed by a statistical method, and some indexes called feature quantities are calculated to detect the abnormality of the device.
  • the sensor data acquisition unit 1001 acquires sensor signals sequentially generated from the sensor 1010 in time series and generates digitized sensor data.
  • the generated sensor data is passed to at least one of the internal variable calculation unit 1003 and the feature quantity calculation unit 1004.
  • the sensor data acquisition unit 1001 may execute a filter process for removing noise irrelevant to the state of the mechanical device 1008, if necessary.
  • the functions and operations of the internal variable holding unit 1002, the internal variable calculation unit 1003, the feature amount calculation unit 1004, and the initialization processing unit 1005 will be described in detail by taking some types of feature amounts as examples.
  • the time when the acquisition of the sensor data used for the feature amount calculation is started is called the "first time point".
  • a total of N sensor data can be obtained from the first time point to the Nth time point.
  • N is an integer of 2 or more.
  • the time point at which the jth sensor data is acquired is referred to as "the jth time point”.
  • j is an integer from 1 to N-1.
  • the sensor data obtained at the j-th time point is referred to as "x j ".
  • the sensor data will be acquired at predetermined time intervals. Needless to say, the sensor data may be acquired irregularly.
  • the internal variable holding unit 1002 holds the internal variables one time before.
  • the number of internal variables is not limited to one, and a plurality of types may be retained depending on the feature amount. The smaller the type of internal variable to be retained is smaller than the number of time series of sensor data on which the feature amount is calculated, the higher the memory reduction effect is.
  • the internal variable holding unit 1002 may hold the feature amount itself as one kind of internal variable.
  • the initialization processing unit 1005 executes the initialization processing for setting the retained internal variables to the initial values.
  • the initial value is a value set at the time of initialization processing.
  • the initialization processing unit 1005 performs the initialization processing at the first time point when the power of the information processing apparatus 1000 is turned on and the sensor data acquisition unit 1001 first acquires the sensor data.
  • the initial value is easily set to 0 as described later, but it does not necessarily have to be 0.
  • an upper limit value and a lower limit value may be set in the vicinity of an appropriate initial value, and the initial value may be appropriately determined between the upper limit value and the lower limit value.
  • the calculation may not be stable if the internal variables and the initial values of the features are set to extremely small values.
  • the calculation is stabilized by setting the initial value to a value other than 0.
  • features whose steady-state values for sensor data are not zero.
  • the convergence of the calculation can be accelerated by setting the value of the assumed feature amount as the initial value.
  • the kurtosis which is one of the feature quantities shown in the present embodiment, has a value of about 3 with respect to the sensor data having a property close to a normal distribution. Therefore, by setting the initial value of kurtosis to 3, it is possible to perform calculations with quick convergence for many types of sensors.
  • the average value mj for the sensor data x1 to xj from the first time point to the jth time point is defined by the following known mathematical formula (1).
  • This formula (1) is also called an arithmetic mean value or an arithmetic mean value.
  • the sensor data at the i -th time point is xj. Therefore, the average value m 3 at the third time point can be expressed by the following mathematical formula (2).
  • the time-series sensor data from the first time point to the third time point is not stored in the memory 1292, and the information of the third time point is held as some internal variables. Then, consider a method of obtaining the average value of the fourth time point based on the internal variables of the third time point and the sensor data of the fourth time point.
  • the average value m 4 at the fourth time point can be expressed as the following mathematical formula (3) according to the mathematical formula (1) which is the definition formula.
  • the sensor data of the time series from the first to the fourth time point is used, and a memory for holding the sensor data for the number of past samples is required. It becomes.
  • the average value m j + 1 at the time of the j + 1 can be expressed by the following mathematical formula (5).
  • the average value m j + 1 at the time point j + 1 is sequentially obtained by the following mathematical formula (6).
  • variable L j + 1 at the time point j + 1 is sequentially obtained by the following mathematical formula (7).
  • Holding the variable at the jth time point in the above formula (7) is equivalent to holding the information of the time j during the sequential calculation.
  • the sensor data x j + 1 is newly acquired at the time point j + 1.
  • the internal variables to be held at the j-th time point in order to obtain the mean value m j + 1 at the j + 1th time point are the mean value m j and the variable L j .
  • the variance v j for the sensor data x 1 to x j from the first time point to the jth time point is defined by the following known mathematical formula (8).
  • an unbiased variance with a denominator of j-1 may be used, but in that case as well, it can be derived by the same procedure.
  • the variance v j + 1 for the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be expressed by the following mathematical formula (11).
  • the variance v j + 1 can be expressed by the following mathematical formula (12).
  • the variable L j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (7). Further, the average value m j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (6). Then, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the j-th time point in order to obtain the variance v j + 1 at the j + 1th time point are the variable L j , the variance v j , and the mean value m j .
  • the standard deviation s j may be held instead of the variance v j .
  • the initial values L 1 and v 1 of the variables L and the variance v should be set to 0 for the sake of simplicity, or the sensor during the initialization process should be used to avoid fluctuations immediately after the initialization process. It is better to decide according to the value of the data. However, if the initial value v 1 of the variance v is set to a value close to 0, the calculation is not stable when the variance v appears in the denominator in the sequential calculation of another feature amount described later. Therefore, it may be good to use a value larger than 0 as the initial value.
  • the standard deviation s j with respect to the sensor data x 1 to x j from the first time point to the jth time point is defined by the following known mathematical formula (15).
  • the standard deviation s j is immediately obtained from the variance v j sequentially obtained by the above procedure. Therefore, the internal variables to be retained for sequentially calculating the standard deviation s j are the same as the variance v j .
  • the root mean square rj for the sensor data x1 to xj from the first time point to the jth time point is defined by the following known mathematical formula (16). This is also called RMS (Root Mean Square: effective value).
  • the root mean square r j + 1 for the sensor data from the first time point to the j + 1st time point can be expressed by the following mathematical formula (17).
  • variable L j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (7). Further, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1. As a result, the internal variables to be held at the j-th time point to obtain the root mean square square root r j 2 at the j + 1th time point are the root mean square square root square r j 2 and the variable L j .
  • the calculation procedure for the sequential calculation of the skewness w is shown.
  • the sensor data x is regarded as a random variable and the mean value m and the standard deviation s of the random variable x are used, the skewness w is defined by the following known mathematical formula (20).
  • the modification of the above formula (21) utilizes the property that the expected value E [x] of the sensor data x is equal to the average value m.
  • variable Aj at the jth time point can be expressed by the following mathematical formula (23).
  • variable A j + 1 at the time point j + 1 can be expressed by the following mathematical formula (24).
  • variable A j + 1 can be expressed as the following formula (25).
  • variable A j + 1 at the time point j + 1 is sequentially obtained by the following formula (26).
  • variable L j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (7). Further, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the j-th time point in order to obtain the variable A j + 1 at the j + 1-th time point are the variable A j and the variable L j .
  • variable B j at the time j can be expressed by the following mathematical formula (27).
  • the skewness w j + 1 with respect to the sensor data x j + 1 from the first time point to the j + 1st time point can be obtained by the following mathematical formula (29).
  • variable A j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (26). Further, the average value m j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (6). Then, the root mean square r j 2 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (19). Further, the standard deviation s j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (15). Then, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the jth time point are the variables A j , Lj, the mean value m j , and the root mean square.
  • the root mean square r j 2 and the standard deviation s j are the variables A j , Lj, the mean value m j , and the root mean square.
  • the root mean square r j 2 and the standard deviation s j are the variables A j , Lj, the mean value m j , and the root mean square.
  • the root mean square r j 2 and the standard deviation s j are the variance v j may be retained instead of the standard deviation s j .
  • the root mean square r j may be held instead of the root mean square r j 2 .
  • the initial values of the variables A j , Lj, the mean value m j , the root mean square square r j 2 , and the standard deviation s j are set to 0. Or, in order to avoid fluctuation immediately after the initialization process, it is better to determine according to the value of the sensor data at the time of the initialization process.
  • the kurtosis k is defined by the following known mathematical formula (30).
  • E [] represents the expected value of the random variable in square brackets.
  • the kurtosis for a random variable having a normal distribution is 3.
  • the modification of the above formula (31) utilizes the property that the expected value E [x] of the sensor data x is equal to the average value m.
  • the expected value E [x 4 ] of the first numerator on the right side is represented by the variable C j .
  • the expected value of the second term and the expected value of the third term of the molecule on the right side can be expressed by variables A j and B j , respectively. Therefore, the kurtosis k j can be expressed by the following mathematical formula (32).
  • the standard deviation s j , the variables A j , and B j can be sequentially calculated by the above formulas (15), formula (26), and formula (27), respectively.
  • variable Cj at the time j can be expressed by the following mathematical formula (33).
  • variable C j + 1 at the time point j + 1 can be expressed by the following mathematical formula (34).
  • variable C j + 1 can be expressed as the following formula (35).
  • the variable B j agrees with the above-mentioned root mean square square square r j 2 . Therefore, the kurtosis k j for the sensor data x 1 to x j from the first time point to the jth time point can be obtained by the following mathematical formula (37).
  • the kurtosis k j + 1 for the sensor data x 1 to x j + 1 from the first time point to the j + 1st time point can be obtained by the following mathematical formula (38).
  • variable C j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (36). Further, the variable A j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (26). Then, the average value m j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (6). Further, the root mean square square r j 2 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (19). Further, the standard deviation s j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (15).
  • the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the jth time point are the variables A j , L j , C j , and the mean value m. j , the root mean square squared r j 2 , and the standard deviation s j .
  • the variance v j may be retained instead of the standard deviation s j .
  • the root mean square r j may be held instead of the root mean square r j 2 .
  • the initial values of the variables A j , L j , C j , the mean value m j , the root mean square square r j 2 , and the standard deviation s j It is better to set to 0 or to determine according to the value of the sensor data at the time of the initialization processing in order to avoid the fluctuation immediately after the initialization processing. Further, it is known that the kurtosis has a value of about 3 when it is known that the sensor data has a property close to a normal distribution. Therefore, the initial value k 1 of the kurtosis k may be set to 3.
  • the maximum value aj for the sensor data x1 to xj from the first time point to the jth time point is defined by the following known mathematical formula (39).
  • the sequential calculation method of the maximum value a can be easily realized by the following method.
  • the maximum value a j at the jth time point is retained, and the larger of the sensor data xj at the j + 1th time point and the maximum value aj at the jth time point is held as the maximum value aj + 1 at the j + 1th time point. .. That is, the maximum value a j + 1 at the time point j + 1 can be sequentially calculated by the following mathematical formula (40).
  • the minimum value n j for the sensor data x 1 to x j from the first time point to the jth time point is defined by the following known mathematical formula (41).
  • the sequential calculation method of the minimum value n can be easily realized by the following method.
  • the minimum value n j at the jth time point is retained, and the smaller of the sensor data xj at the j + 1th time point and the minimum value nj at the jth time point is held as the minimum value nj + 1 at the j + 1th time point. .. That is, the minimum value n j + 1 at the time point j + 1 can be sequentially calculated by the following mathematical formula (42).
  • the peak value pj from the first time point to the jth time point is defined by the following formula (43).
  • the peak value is the largest absolute value among the sensor data x1 to xj .
  • the sequential calculation method of the peak value can be realized by the following method.
  • the peak value p j at the jth time point is retained, and the larger of the sensor data x j at the j + 1th time point and the peak value pj at the jth time point is held as the peak value pj + 1 at the j + 1th time point. That is, the peak value p j + 1 at the time point j + 1 can be sequentially calculated by the following mathematical formula (44).
  • the peak value is obtained by the following formula (45). You can also.
  • the peak value can be obtained by the following mathematical formula (46).
  • the internal variable to be held at the jth point in order to obtain the peak value pj + 1 at the j + 1th time point is the peak value pj .
  • the maximum value a j and the minimum value n j may be held instead of the peak value p j .
  • Peak The peak value is the difference between the maximum value and the minimum value. It is also called peak-to-peak.
  • the peak peak value pp j for the sensor data x1 to xj from the first time point to the jth time point is defined by the following known mathematical formula (47).
  • the peak peak value pp j can be immediately obtained by sequentially calculating the maximum value a j and the minimum value n j from the above-mentioned mathematical formulas (40) and (42).
  • the peak peak value pp j + 1 can be obtained by the following known mathematical formula (48).
  • the internal variables to be held at the jth time point in order to obtain the peak peak value pp j + 1 at the j + 1th time point are the maximum value a j and the minimum value n j .
  • the wave height rate pr j with respect to the sensor data x1 to xj from the first time point to the jth time point is defined by the following mathematical formula (49).
  • the peak factor pr j is a value obtained by dividing the above-mentioned peak value p j by the above-mentioned root mean square r j .
  • the crest factor pr j is sometimes referred to as peak to RMS (peak-to-RMS) or crest factor. If it is known that the sensor data has a value biased in the positive direction, the maximum value a j may be used instead of the peak value p j .
  • the wave height rate pr j + 1 for the sensor data x1 to xj from the first time point to the j + 1st time point can be obtained by the following mathematical formula (50).
  • the peak value pj + 1 is sequentially calculated using the above-mentioned formula (44) or formula (45), and the root mean square r j + 1 is sequentially calculated using the formula (19).
  • the wave height rate pr j + 1 can be obtained immediately.
  • the internal variables to be held at the jth time point are the peak value pj and the root mean square square r j2 .
  • the maximum value a j and the minimum value n j may be held instead of the peak value p j .
  • the root mean square r j 2 of the root mean square instead of the root mean square r j may be retained.
  • the internal variable calculation unit 1003 sequentially calculates the internal variables for calculating each feature amount.
  • the feature amount calculation unit 1004 can determine the feature amount corresponding to the sensor data x1 to xj from the first time point to the Nth time point.
  • the updated formula and the internal variables held by the internal variable holding unit 1002 are different for each type of feature amount.
  • the feature amount is that the internal variable holding unit 1002 holds the internal variable, the internal variable calculation unit 1003 sequentially calculates the internal variable, and the feature amount calculation unit 1004 calculates the feature amount based on the internal variable. It is the same regardless of the type of.
  • FIG. 6 is a diagram showing time-series waveforms of the motor torque and various feature quantities in the first embodiment.
  • the notation of the symbol attached to the feature amount and the like will be omitted as appropriate.
  • FIG. 6 (a) shows the same motor torque waveform as that shown in FIG. 5 (b).
  • the horizontal axis of FIGS. 6A and 6B represents time.
  • FIG. 6A shows a time series of mean value, standard deviation, root mean square (RMS), maximum value and minimum value having the same unit [Nm] as the motor torque among the dimensional features.
  • the waveform is shown.
  • the peak value is the difference between the maximum value and the minimum value, and is not shown for simplicity. Further, the peak value is a value indicating a value having a large absolute value among the maximum value and the minimum value, and is not shown for simplicity.
  • the variance is the square of the standard deviation and is not shown because the units are different.
  • FIG. 6B shows a time-series waveform of skewness, kurtosis, and crest factor, which are dimensionless features, that is, features having no unit.
  • the average value was almost 0 in the period Ts1, gradually increased in the period Ta, slightly decreased in the period Te, gradually decreased in the period Td, and slightly decreased in the period Ts2.
  • the standard deviation is a small value that is almost constant in the period Ts1, gradually increases in the period Ta, slightly decreases in the period Te, gradually increases in the period Td, and slightly decreases in the period Ts2.
  • the root mean square is a small value that is almost constant in the period Ts1, gradually increases in the period Ta, decreases slightly in the period Te, is almost constant in the period Td, and decreases slightly in the period Ts2.
  • the maximum value is a small positive value that is almost constant in the period Ts1, increases every time the motor torque exceeds the past maximum value in the period Ta, and is constant in the period Te, the period Td, and the period Ts2.
  • the minimum value is a small negative value that is almost constant in the period Ts1, constant in the period Ta and Te, decreases every time the motor torque exceeds the past minimum value in the period Td, and is constant in the period Ts2. ..
  • the skewness is almost 0 in the period Ts1, increases sharply to a large positive value at the time Tr1 at which the period Ta starts, gradually decreases to a small value in the subsequent period Ta, and increases slightly in the period Te. It decreased slightly in the period Td and increased slightly in the period Ts2.
  • the kurtosis increases from about 2 to about 3 in the period Ts1, increases sharply to a large positive value at the time Tr1 when the period Ta starts, gradually decreases to a small value in the subsequent period Ta, and slightly in the period Te. In the period Td, it increased slightly and then decreased, and in the period Ts2, it increased slightly.
  • the wave height rate increases from about 2 to about 3 in the period Ts1, increases sharply to a large positive value at the time Tr1 when the period Ta starts, gradually decreases to a small value in the subsequent period Ta, and slightly in the period Te. In the period Td, it is almost constant, and in the period Ts2, it increases slightly.
  • the feature amount at each time point shown in FIGS. 6A and 6B is a sequential calculation of the feature amount with respect to the motor torque in the past from each time point. Therefore, for example, the average value is a value obtained by batch calculation while holding the motor torque for the entire period (Ts1 + Ta + Te + Td + Ts2) in FIG. 6 (a), and a time Tr5 which is the final time in FIG. 6 (a) sequentially obtained. Matches the value of.
  • the feature amount calculation unit 1004 in order to make the transition of the feature amount easy to change, the feature amount calculation unit 1004 also calculates each feature amount for each sensor data acquisition cycle.
  • the internal variable calculation unit 1003 needs to sequentially calculate the internal variable for each acquisition cycle of the sensor data.
  • the feature amount calculation unit 1004 is necessary. The feature amount may be calculated accordingly.
  • the feature amount calculation unit 1004 may be configured to calculate the feature amount only once at the time point Tr5, that is, when the data score N of the sensor data reaches 12000. With this configuration, it is possible to reduce the computational load of the information processing apparatus 1000.
  • FIG. 7 is a diagram showing an example of the secular variation of the feature amount in the first embodiment.
  • FIG. 7 shows the results of periodically plotting the features calculated by the above method. Regular means, for example, monthly.
  • the horizontal axis of FIG. 7 is time, which represents the operating time of the mechanical device 1008.
  • the state diagnosis unit 1006 sets a threshold value Fth1 for the feature amount, and when the feature amount exceeds the threshold value Fth1, diagnoses that an abnormality has occurred in the mechanical device 1008.
  • the feature amount any of the above-mentioned feature amounts may be used. Further, a plurality of feature quantities may be calculated and individual threshold values may be set for each feature quantity.
  • the time Tta0 is the time when the operation is started.
  • the feature amount from the time Tta0 to the time before the time Tta1 keeps a value smaller than the threshold value Fth1 although there is some variation.
  • the state diagnosis unit 1006 makes a diagnosis that the mechanical device 1008 is normal, and outputs the diagnosis result.
  • the state diagnosis unit 1006 makes a diagnosis that an abnormality has occurred in the mechanical device 1008 at the time Tta1, and outputs the diagnosis result.
  • the threshold value Fth1 may be determined based on the feature amount when an abnormality occurs in another mechanical device in the past. Further, the threshold value Fth1 may be determined with reference to the feature amount immediately after the operation of the mechanical device 1008 is started. Further, the threshold value Fth1 may be set dynamically. For example, while operating a plurality of mechanical devices 1008, the feature amounts thereof may be calculated individually, and the threshold value Fth1 may be changed periodically in consideration of the variation in the feature amounts for each device. Further, the threshold value Fth1 may be determined based on the feature amount obtained by a simulation or the like in consideration of the characteristics of the mechanical device 1008.
  • the threshold value Fth1 represents the upper limit value of the normal state, but when a feature amount in which a smaller value indicates an abnormality is used, the lower limit value of the normal state is used as the threshold value. It may be set as Fth1. Further, two threshold values consisting of an upper limit value and a lower limit value may be used in combination. Further, the diagnosis result does not have to be one, and the diagnosis may be performed for a plurality of feature quantities and the diagnosis result may be output for each feature quantity.
  • the method is not limited to this method.
  • a statistical method such as performing a test based on the distribution of the degree of abnormality for a certain period may be adopted.
  • FIG. 8 is a flowchart provided for explaining the information processing method according to the first embodiment.
  • the sensor data acquisition unit 1001 acquires the measured value of the physical quantity as sensor data from the sensor 1010 that measures the physical quantity of the mechanical device 1008 at each time point from the first time point to the Nth time point (N is an integer of 2 or more).
  • the internal variable holding unit 1002 holds a number of internal variables smaller than N, which are sequentially calculated based on the sensor data (step S102).
  • the internal variable calculation unit 1003 calculates the internal variable corresponding to the j + 1th time point (j is an integer from 1 to N-1) based on the sensor data at the j + 1st time point and the internal variable corresponding to the jth time point. (Step S103).
  • the feature amount calculation unit 1004 calculates the feature amount obtained by extracting the statistical features included in the sensor data from the first time point to the Nth time point based on the internal variables of the Nth time point (step S104).
  • the state diagnosis unit 1006 diagnoses the state of the mechanical device based on the feature amount (step S105).
  • a step of executing an initialization process for determining an internal variable at a first time point to a value between a preset maximum value and a minimum value may be included.
  • the sensor data acquisition unit acquires the measured value of the physical quantity of the mechanical device measured by the sensor, and starts from the first time point among the measured values.
  • the measured values up to the Nth time point (N is an integer of 2 or more) are retained as sensor data.
  • the internal variable holding unit holds a number of internal variables smaller than N, which are sequentially calculated in time series based on the sensor data.
  • the internal variable calculation unit calculates the internal variable corresponding to the j + 1th time point (j is an integer from 1 to N-1) based on the sensor data at the j + 1st time point and the internal variable corresponding to the jth time point.
  • the feature amount calculation unit calculates the feature amount obtained by extracting the statistical features included in the sensor data from the first time point to the Nth time point based on the internal variables of the Nth time point. According to the information processing apparatus configured as described above, it is not necessary to hold the sensor data from the first time point to the Nth time point by holding the number of internal variables smaller than N. As a result, it is possible to obtain an unprecedented remarkable effect that the feature amount can be calculated without using a computer having a high computing power and a computer having a large storage capacity.
  • the information processing device includes a state diagnosis unit that diagnoses the state of the mechanical device based on the feature amount.
  • a state diagnosis unit that diagnoses the state of the mechanical device based on the feature amount.
  • the information processing apparatus includes an initialization processing unit that executes an initialization process for determining an internal variable at a first time point to a value between a preset maximum value and a preset minimum value. May be.
  • the physical quantity at each time point from the first time point to the Nth time point (N is an integer of 2 or more) from the sensor for measuring the physical quantity of the mechanical device.
  • the measured value of is acquired as sensor data.
  • a smaller number of internal variables than N which are sequentially calculated in time series based on the sensor data, are held.
  • the internal variable corresponding to the j + 1th time point (j is an integer from 1 to N-1) is calculated based on the sensor data at the j + 1st time point and the internal variable corresponding to the jth time point.
  • the feature amount obtained by extracting the statistical features included in the sensor data from the first time point to the Nth time point is calculated based on the internal variables of the Nth time point.
  • the state of the mechanical device is diagnosed based on the feature amount. According to the information processing method including the processes of the first to fifth steps, a smaller number of internal variables than N are held, so that it is necessary to hold the sensor data from the first time point to the Nth time point. There is no. As a result, it is possible to obtain an unprecedented remarkable effect that information processing for diagnosing the state of a mechanical device can be performed without using a computer having a high computing power and a computer having a large storage capacity.
  • the internal variable at the first time point is initially determined to be a value between a preset maximum value and a preset minimum value. It may include a step to execute the conversion process. By including such an initialization processing step, the probability of erroneous determination can be reduced and the accuracy of the diagnostic processing can be improved.
  • FIG. 9 is a diagram showing a configuration example of the information processing system 100A including the information processing apparatus 2000 according to the second embodiment.
  • the information processing device 1000 is replaced with the information processing device 2000
  • the device control unit 1099 is replaced with the device control unit 2099.
  • the sensor data acquisition unit 1001 is replaced by the sensor data acquisition unit 2001
  • the internal variable calculation unit 1003 is replaced by the internal variable calculation unit 2003
  • the feature quantity calculation unit 1004 is replaced by the feature quantity calculation unit 2004. It has been replaced.
  • Other configurations are the same as or equivalent to the information processing system 100 shown in FIG.
  • the same or equivalent components are designated by the same reference numerals, and duplicate description will be omitted.
  • the sensor data acquisition unit 2001 In addition to the processing of the sensor data acquisition unit 1001, the sensor data acquisition unit 2001 generates a calculation permission flag for determining whether or not to update the internal variable based on the operation signal. More generally, the sensor data acquisition unit 2001 has at least one of an acquisition cycle, which is an interval of time for acquiring sensor data, an operation signal indicating the operation status of the mechanical device 1008, and a data value in the sensor data. Two or more time points between the first time point and the jth time point (j is an integer from 1 to N-1 (where N is an integer of 2 or more)) are determined based on the above, and two or more points are determined. Generates a compute permission flag to update internal variables based on the time point in.
  • the internal variable calculation unit 2003 updates the internal variable when the calculation permission flag contains the content of updating the internal variable, in addition to the processing of the internal variable calculation unit 1003. On the other hand, when the calculation permission flag does not include the content to update the internal variable, the internal variable calculation unit 2003 performs a process of inheriting the value of the internal variable from the previous time. More generally, the internal variable calculation unit 2003 uses the internal variable calculation unit 2003 to indicate each u-th time point indicated by the calculation permission flag (u is an integer of 1 or more and j or less (where j is an integer from 1 to N-1 and N). Is an integer of 2 or more)), the internal variable at the u + 1 time point is determined to be the value of the internal variable at the uth time point.
  • the feature amount calculation unit 2004 updates the feature amount when the calculation permission flag includes the content to update the feature amount, in addition to the processing of the feature amount calculation unit 1004. If the calculation permission flag does not include the content to update the feature amount, the feature amount calculation unit 2004 performs a process of inheriting the feature amount from the previous time.
  • the device control unit 2099 outside the information processing device 2000 outputs an operation signal indicating the operation status of the mechanical device 1008 in addition to the processing of the device control unit 1099.
  • the mechanical device 1008 is, for example, a mechanical device driven by a motor
  • sensor data with less noise can be obtained when the feature amount when the speed is constant is used, and the state of the mechanical device 1008 can be estimated with high accuracy. be. Therefore, by configuring the internal variables and the values of the feature amount to be updated only when the speed of the motor 1009 is constant, it becomes possible to calculate the feature amount more useful for the state diagnosis of the mechanical device 1008. ..
  • FIG. 10 is a diagram showing time-series waveforms of motor torque and various feature quantities in the second embodiment.
  • FIG. 10 (a) shows the same waveform as the time-series waveform of the motor speed shown in FIG. 5 (a). Further, FIG. 10B shows a time-series waveform of the calculation permission flag determined by the sensor data acquisition unit 2001.
  • the flag value is set to "1" when the internal variable is updated, and the flag value is set to "0" when the internal variable is not updated.
  • the sensor data acquisition unit 2001 sets the flag value of the calculation permission flag to "1” in the period Te from the time Tr2 to the time Tr3 so that the internal variable is updated when the speed is not "0" and is constant. It is supposed to be. Further, in a period other than the period Te, that is, in the period Ts1, the period Ta, the period Td, and the period Ts2, the sensor data acquisition unit 2001 sets the flag value of the calculation permission flag to “0”. In FIG. 9, an example of constant speed is set to 500 [r / min].
  • the internal variable is not updated when the speed is 0 is that the mechanical device 1008 does not operate when the motor 1009 is stopped, and the information corresponding to the state of the mechanical device 1008 does not appear in the sensor data. Can be mentioned. However, depending on the configuration of the mechanical device 1008, information according to the state of the mechanical device 1008 may appear in the sensor data even when the motor 1009 is stopped. In such a case, the internal variables may be configured to be updated even when the system is stopped. The operation when the flag value of the calculation permission flag is "0" does not have to change the internal variables and the feature amount. For example, even if the calculation is performed so that the current value is overwritten with the previous value. The calculation itself may be skipped.
  • the speed and acceleration of the motor are used as operation signals.
  • the flag value of the calculation permission flag is set to "1".
  • the reference acceleration and speed may be determined in advance with reference to the specifications of the mechanical device 1008 or the motor 1009.
  • the operation signal is obtained from the device control unit 2099.
  • the flag value of the calculation permission flag is set to "1" when the predetermined time is reached after the motor 1009 starts to move, and the flag value of the calculation permission flag is set to "0" after another predetermined time has elapsed. You may try to do it.
  • the sampling cycle which is the time interval for acquiring the sensor data.
  • FIG. 10 (c) shows a time-series waveform of the motor torque similar to that shown in FIG. 5 (b). Further, FIG. 10 (c) shows a time-series waveform having the same feature amount as that shown in FIG. 6 (a). Further, FIG. 10 (d) shows a time-series waveform having the same feature amount as that shown in FIG. 6 (b).
  • the values of the feature amounts in FIGS. 10 (c) and 10 (d) did not change in the period Ts1, the period Ta, the period Td, and the period Ts2 excluding the period Te.
  • the internal variable calculation unit 2003 and the feature amount calculation unit 2004 are configured so that the internal variables and the feature amount are not updated during the period other than the period Te. Because it is.
  • the period Te all the features are updated at any time, and by the time Tr3, which is the end point of the period Te, the features have converged to almost constant values. According to such a configuration, in addition to the effect of the first embodiment, the state of the mechanical device 1008 can be estimated with higher accuracy.
  • the sensor data acquisition unit has an acquisition cycle which is an interval of time for acquiring sensor data, an operation signal indicating the operation status of the mechanical device, and a sensor. Two or more time points between the first time point and the jth time point are determined based on at least one of the data values in the data, and the calculation permission flag is generated based on the determined two or more point points.
  • the internal variable calculation unit determines the internal variable at the u + 1 time point as the value of the internal variable at the uth time point for each uth time point (u is an integer of 1 or more and j or less) indicated by the calculation permission flag. According to the information processing device configured in this way, sensor data with less noise can be obtained according to the operating state of the mechanical device. As a result, after obtaining the effect of the first embodiment, the effect that the state of the mechanical device can be estimated with high accuracy can be obtained.
  • FIG. 11 is a diagram showing a configuration example of the information processing system 100B including the information processing apparatus 3000 according to the third embodiment.
  • the information processing apparatus 2000 is replaced with the information processing apparatus 3000 in the configuration of the information processing system 100A shown in FIG.
  • the sensor data acquisition unit 2001 is replaced by the sensor data acquisition unit 3001
  • the initialization processing unit 1005 is replaced by the initialization processing unit 3005.
  • Other configurations are the same as or equivalent to the information processing system 100A shown in FIG.
  • the same or equivalent components are designated by the same reference numerals, and duplicate description will be omitted.
  • the sensor data acquisition unit 3001 In addition to the processing of the sensor data acquisition unit 1001 shown in FIG. 1, the sensor data acquisition unit 3001 generates an initialization trigger for determining a time point for initializing an internal variable. More generalized and more specifically, the sensor data acquisition unit 3001 has a acquisition cycle which is an interval of time for acquiring sensor data, an operation signal indicating the operation status of the mechanical device 1008, and sensor data. Determine one or more time points between the first and Nth time points based on at least one, and generate an initialization trigger to initialize the internal variables based on the determined one or more time points. ..
  • the initialization processing unit 3005 performs the initialization processing of the internal variables based on the time point indicated by the initialization trigger in addition to the processing of the initialization processing unit 1005.
  • the tendency of the feature amount may differ depending on the operating conditions such as acceleration, rotation at a constant speed, deceleration, and stop. In such a case, it may be possible to estimate the state of the mechanical device 1008 with high accuracy by separating the calculation of the feature amount for each operating condition. Therefore, by initializing the internal variables at the time when the operating condition changes, such as when the speed change of the motor 1009 is changed, it becomes possible to calculate the feature amount which is more useful for estimating the state of the mechanical device 1008.
  • FIG. 12 is a diagram showing time-series waveforms of the motor torque and various feature quantities in the third embodiment.
  • FIG. 12 (a) shows the same waveform as the time-series waveform of the motor speed shown in FIG. 5 (a). Further, FIG. 12B shows a time-series waveform of the initialization trigger determined by the sensor data acquisition unit 3001.
  • the signal level is set to "1" when the internal variable is initialized, and the signal level is set to "0" when the internal variable is not initialized.
  • the sensor data acquisition unit 3001 sets the signal level of the initialization trigger at time Tr1, time Tr2, time Tr3 and time Tr4 to "1" so as to initialize the internal variable when the motor speed changes, that is, when acceleration occurs.
  • the signal level is set to "0" at other times.
  • the jerk of the motor that is, the time derivative of the acceleration is used as the operation signal.
  • the signal level of the initialization trigger is set to "1".
  • the reference jerk size should be determined in advance with reference to the specifications of the mechanical device 1008 or the motor 1009.
  • the operation signal is obtained from the device control unit 2099.
  • the signal level of the initialization trigger may be momentarily set to "1" when the predetermined time is reached after the motor 1009 starts to move.
  • the signal level of the initialization trigger may be momentarily set to "1" every time a specific time elapses.
  • the sampling cycle which is the time interval for acquiring the sensor data.
  • the initialization trigger may be determined so that the initialization process is performed when the motor 1009 starts to move.
  • FIG. 12 (c) shows a time-series waveform of the motor torque similar to that shown in FIG. 5 (b). Further, FIG. 12 (c) shows a time-series waveform having the same feature amount as that shown in FIG. 6 (a). Further, FIG. 12 (d) shows a time-series waveform having the same feature amount as that shown in FIG. 6 (b).
  • each feature quantity of FIGS. 10 (c) and 10 (d) begins to change immediately after the time Tr0, the time Tr1, the time Tr2, the time Tr3 and the time Tr4, and each period Ts1, period Ta, period Te, and period By the end of Td and the period Ts2, each feature quantity has converged to a substantially constant value.
  • the state of the mechanical device 1008 can be estimated with higher accuracy.
  • it can be applied to the case where the motor 1009 is operated under the condition that the motor 1009 does not have a constant speed, so that the applicable range can be expanded.
  • the sensor data acquisition unit has an acquisition cycle which is an interval of time for acquiring sensor data, an operation signal indicating the operation status of the mechanical device, and a sensor.
  • an acquisition cycle which is an interval of time for acquiring sensor data, an operation signal indicating the operation status of the mechanical device, and a sensor.
  • the initialization processing unit executes a process of initializing an internal variable based on an initialization trigger. According to the information processing apparatus configured in this way, it is possible to separate the calculation of the feature amount for each operating condition.
  • the effect of the first embodiment it is possible to further obtain the effect that the state of the mechanical device can be estimated with high accuracy. Further, since it can be applied to the case where the motor is operated under the condition that the motor does not have a constant speed, the effect that the application range regarding the operating condition can be further expanded after obtaining the effect of the second embodiment can be obtained.
  • FIG. 13 is a diagram showing a configuration example of the information processing system 100C including the information processing apparatus 4000 according to the fourth embodiment.
  • the information processing apparatus 3000 is replaced with the information processing apparatus 4000 in the configuration of the information processing system 100B shown in FIG.
  • the internal variable calculation unit 1003 is replaced with the internal variable calculation unit 4003
  • the feature amount calculation unit 1004 is replaced with the feature amount calculation unit 4004
  • the state diagnosis unit 1006 is replaced with the state diagnosis unit 4006.
  • a convergence degree calculation unit 4007 is newly provided.
  • Other configurations are the same as or equivalent to the information processing system 100B shown in FIG.
  • the same or equivalent components are designated by the same reference numerals, and duplicate description will be omitted.
  • the internal variable calculation unit 4003 calculates the internal variable based on the oblivion coefficient in addition to the processing of the internal variable calculation unit 1003.
  • the feature amount calculation unit 4004 also calculates the feature amount based on the forgetting coefficient in addition to the processing of the feature amount calculation unit 1004.
  • the forgetting coefficient is a coefficient for weighting the feature amount so that the effect of the sensor data with the new time on the feature amount is larger than the effect of the sensor data with the old time on the feature amount. It has a small value.
  • the convergence degree calculation unit 4007 calculates the degree of convergence, which is an index that quantitatively indicates the degree of convergence in the calculation of the feature amount, based on the internal variables.
  • the state diagnosis unit 4006 diagnoses the state of the mechanical device 1008 based on the feature amount at the time when the convergence degree calculated by the convergence degree calculation unit 4007 satisfies a certain condition.
  • the calculation procedure for the sequential calculation of the exponential moving average is shown.
  • a weighted average value considering the weight given for each time is known.
  • the weighted average value there is a method of averaging the weighted sensor data weighted by multiplying the forgetting coefficient ⁇ .
  • the forgetting coefficient ⁇ is a coefficient for weighting the time so that the effect of the time on the new sensor data is greater than the effect of the time on the old sensor data.
  • this process is referred to as an exponential moving average, and the value is referred to as an exponential moving average value.
  • the exponential moving average is sometimes referred to as the exponential smoothing moving average.
  • "1- ⁇ " obtained by subtracting the forgetting coefficient ⁇ from "1" may be referred to as a smoothing coefficient.
  • the exponential moving average value m'j for the sensor data x1 to xj from the first time point to the jth time point is defined by the following mathematical formula (51).
  • the exponential moving average value m'3 at the third time point can be expressed by the following mathematical formula (52).
  • the coefficient in the latest data is larger than that in the old data as described above (0.81 ⁇ 0.9 ⁇ 1). It can be said that the recent data, that is, the new data of the time, is emphasized when calculating the average.
  • the exponential moving average value m'j + 1 with respect to the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be transformed by the following mathematical formula (53).
  • the exponential moving average value m'j + 1 at the time point j + 1 is sequentially obtained by the following mathematical formula (54).
  • the sensor data x j + 1 is newly acquired at the time point j + 1. Further, the variable L' j at the time of the j + 1 can be expressed by the following mathematical formula (55).
  • variable L' j + 1 at the time of the j + 1 can be expressed by the following mathematical formula (56).
  • variable L'j + 1 at the time point j + 1 can be sequentially calculated by the following mathematical formula (57).
  • the internal variables that should be retained for calculating the exponential moving average m'j + 1 at the j + 1th time point are the exponential moving average value m'j and the variable L' j at the jth time point.
  • the forgetting coefficient ⁇ is preferably set to a value larger than 0 and smaller than 1.
  • the forgetting coefficient ⁇ is close to 1, it is difficult to forget the past information, and when it is close to 0, it is easy to forget the past information.
  • the forgetting coefficient ⁇ is set to 1, the information of the sensor data is not forgotten, and the formula of the sequential calculation is the same as that of the first embodiment. Further, when the forgetting coefficient is set to 0, the information of all the sensor data is forgotten when one time elapses, so that the feature amount for the sensor data of the past two or earlier cannot be obtained.
  • the forgetting factor should be set to a value close to 1 and less than 1, for example, 0.9 to 0.999.
  • the time constant ⁇ can be calculated by the following mathematical formula (58).
  • the oblivion coefficient ⁇ may be set by using the estimated frequency (unit rad / s) which is the reciprocal of the time constant ⁇ .
  • the calculation procedure for the sequential calculation of the exponential movement variance is shown.
  • the variance using the weighted sensor data weighted by multiplying the forgetting coefficient ⁇ for each time is called “exponential movement variance”.
  • the exponential movement variance v'j + 1 for the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be expressed by the following mathematical formula (60).
  • variable L'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (57). Further, the exponential moving average value m'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (54). Then, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the j-th time point in order to obtain the exponential moving variance v'j + 1 at the j + 1th time point are the variable L' j , the exponential moving variance v'j , and the exponential moving average value m'j .
  • the exponential movement standard deviation s'j described later may be retained.
  • the initial values L' 1 , v'1, m'1 of the variables L', exponential moving variance v'and exponential moving average m' are set to 0 .
  • the exponential movement standard deviation s'j with respect to the sensor data x1 to xj from the first time point to the jth time point is defined by the following formula (64) from the exponential movement variance v'j .
  • the exponential movement standard deviation s'j is immediately obtained from the exponential movement variance v'j sequentially obtained by the above procedure. Therefore, the internal variable to be retained for sequentially calculating the exponential movement standard deviation s'j is the same as the exponential movement variance v'j .
  • the calculation procedure for the sequential calculation of the exponentially moving root mean square is shown.
  • the root mean square using the weighted sensor data weighted by multiplying the forgetting coefficient ⁇ for each time is called an "exponentially moving root mean square".
  • the exponential movement root mean square r'j for the sensor data x1 to xj from the first time point to the jth time point is defined by the following mathematical formula (65).
  • the exponential root mean square r'j + 1 for the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be expressed by the following mathematical formula (66).
  • variable L'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (57). Further, the sensor data x j + 1 is newly acquired at the time point j + 1. As a result, the internal variables to be held at point j in order to obtain the square root of the exponentially moved root mean square at the j + 1 point are the square r'j 2 of the exponentially moved root mean square and the variable L' j . be.
  • the calculation procedure for the sequential calculation of the exponential movement skewness is shown.
  • the skewness using the weighted sensor data weighted by multiplying the forgetting coefficient ⁇ for each time is called "exponential skewness”.
  • the exponential movement skewness w'j for the sensor data x1 to xj from the first time point to the jth time point is the following formula (69). Can be expressed by.
  • the exponential moving average value m'j of the above formula (69) can be sequentially calculated by the above formula (54), and the exponential moving standard deviation s'j of the above formula (69) can be calculated sequentially by the above formula (54). It can be calculated sequentially in (63) and (64).
  • variables A'j and B'j at the jth time point can be expressed by the following mathematical formulas (70) and (71).
  • variable A'j + 1 at the time of the j + 1 can be expressed by the following mathematical formula (72).
  • variable A'j + 1 can be expressed as the following formula (73).
  • variable A'j + 1 at the time j + 1 is sequentially obtained by the following formula (74).
  • the internal variables to be held at the jth time point in order to obtain the variable at the j + 1th time point are the variables L' j and A'j .
  • the variable L'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (57).
  • the variable B'j at the jth time point can be expressed by the following mathematical formula (75).
  • variable B'j matches the square root r'j 2 of the exponential moving square root mean square as in the above formula (75), which is sequential in the above formula (68). Is required. Therefore, the exponential movement skewness w'j for the sensor data x1 to xj from the first time point to the jth time point can be obtained by the following mathematical formula (76).
  • the exponential movement skewness w'j + 1 with respect to the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be obtained by the following mathematical formula (77).
  • variable A'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (74). Further, the exponential moving average value m'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (54). Then, the square root r'j 2 of the exponential moving square root mean square at the j + 1st time point is sequentially obtained by the above-mentioned mathematical formula (68). Further, the exponential movement standard deviation s'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (64). Then, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1.
  • the internal variables to be held at the jth time point are the variables A'j, L' j , and the exponent.
  • the exponential movement variance v'j may be retained instead of the exponential movement standard deviation s'j .
  • the exponential moving square root mean square r'j may be held instead of the square root r'j 2 of the exponentially moving root mean square.
  • the calculation procedure for the sequential calculation of the exponential movement kurtosis is shown.
  • the kurtosis using the weighted sensor data weighted by multiplying the forgetting coefficient ⁇ for each time is called "exponential kurtosis”.
  • the exponential movement kurtosis k'j for the sensor data x1 to xj from the first time point to the jth time point is the following formula (78). Can be expressed by.
  • variable C'j at the jth time point can be expressed by the following mathematical formula (79).
  • variable C'j + 1 at the time of the j + 1 can be expressed by the following mathematical formula (80).
  • variable C'j + 1 can be expressed as the following formula (81).
  • the variable B'j matches the square root r'j 2 of the exponentially moving root mean square, and the square r'j 2 of the exponentially moved root mean square is described above. It is sequentially obtained by the formula (68) of. Therefore, the exponential movement kurtosis k'j for the sensor data x1 to xj from the first time point to the jth time point can be obtained by the following mathematical formula (83).
  • the exponential movement kurtosis k'j + 1 with respect to the sensor data x1 to xj + 1 from the first time point to the j + 1st time point can be obtained by the following mathematical formula (84).
  • variable C'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (82). Further, the variable A'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (74). Then, the exponential moving average value m'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formula (54). Further, the square root r'j + 1 2 of the exponential moving square root mean square at the j + 1st time point is sequentially obtained by the above-mentioned mathematical formula (68).
  • the exponential movement standard deviation s'j + 1 at the time point j + 1 is sequentially obtained by the above-mentioned mathematical formulas (63) and (64). Then, the sensor data x j + 1 at the time j + 1 is newly acquired at the time j + 1. As a result, the square r'j 2 of the first time point or the exponential moving square root mean square is sequentially obtained by the above formula (68).
  • the internal variables to be held at the jth time point are the variables A'j, L' j , and the exponential movement average value m'j .
  • Exponential movement square root mean square r'j 2 Exponential movement square root mean square r'j 2 and exponential movement standard deviation s'j .
  • the exponential movement variance v'j may be retained instead of the exponential movement standard deviation s'j .
  • the exponential moving square root mean square r'j may be held instead of the square root r'j 2 of the exponentially moving root mean square.
  • the exponential movement kurtosis k'j has a value of about 3 when it is known that the sensor data has a property close to a normal distribution. Therefore, the initial value k ' 1 of the exponential movement kurtosis k'j may be set to 3.
  • the calculation of the maximum value, the minimum value, the peak peak value, the peak value, and the peak rate (hereinafter, appropriately referred to as "maximum value, etc.") in which the latest time is weighted is as shown in the first embodiment. It is not possible to update the current maximum value, etc. from the previous maximum value, etc. However, for example, some candidate values that are likely to be the maximum value may be retained as internal variables, and the retained candidate values may be deleted after the specified time has elapsed. By doing so, it is possible to calculate the maximum value or the like corresponding to the sensor data of the latest time.
  • the degree of convergence that quantitatively represents the degree of convergence of the sequential calculation can be calculated from the internal variables.
  • the variable L' j which is one of the internal variables, is sequentially calculated by the mathematical formula (57), it asymptotically approaches the value of 1 / (1- ⁇ ) regardless of the input sensor data. You can see that. Therefore, when the initial value L' 1 of the variable L' j is 0 or more and less than 1 / (1- ⁇ ), the following formula (85) can be used so that the degree of convergence can be a value from 0 to 1. It is conceivable to define the degree of convergence di with.
  • the degree of convergence is calculated using the above formula (85), but the degree of convergence is not limited to the above, and any definition can be used as long as the value increases with the convergence of the calculation. good. However, the value increases monotonically as the calculation converges, and the one that asymptotics to a specific value is easy to use as a reference.
  • ⁇ (1- ⁇ ) L' j ⁇ M which is obtained by multiplying the entire right side of the above formula (85) by the Mth power (M is a positive real number), is defined as the degree of convergence.
  • This degree of convergence is a preferable definition expression because it takes a value from 0 to 1 and increases monotonically with the convergence of the calculation.
  • FIG. 14 is a diagram showing time-series waveforms of motor torque and various feature quantities in the fourth embodiment.
  • FIG. 14 (a) shows the same waveform as the time-series waveform of the motor speed shown in FIG. 5 (a). Further, FIG. 14 (b) shows the same waveform as the time-series waveform of the initialization trigger shown in FIG. 12 (b).
  • FIG. 14C shows a time-series waveform of the degree of convergence calculated by the degree of convergence calculation unit 4007 using the mathematical formula (85).
  • the degree of convergence is 0 at the time Tr0 when the sequential calculation starts, the time Tr1, the time Tr2, the time Tr3, and the time Tr4 where the initialization trigger indicates “1”. This is because the initialization process is executed and the variable L' j at the jth time point, which is one of the internal variables, becomes 0.
  • the j-th time point indicates the time point when the initialization trigger indicates "1".
  • the degree of convergence has increased since each initialization process was performed, and gradually approaches 1.
  • FIG. 14 (d) shows a time-series waveform of the motor torque similar to that shown in FIG. 5 (b). Further, FIG. 14D shows a time-series waveform of the exponential moving average value, the exponential moving standard deviation, and the exponential moving square root mean square (RMS) having the same unit [Nm] as the motor torque. The notation of "exponential movement" is omitted to avoid complicating the figure.
  • FIG. 14 (e) shows a time-series waveform of the exponential movement skewness and the exponential movement kurtosis, which are dimensionless features, that is, features having no unit.
  • the forgetting coefficient ⁇ is 0.975.
  • FIG. 14D the notation of “exponential movement” is omitted.
  • the exponential moving average value follows the motor torque over the entire period from time Tr0 to time Tr5, although it is slightly behind the motor torque. This shows the nature of the present embodiment in which the average value is calculated by giving a large weight to the latest short-time motor torque.
  • the exponential movement standard deviation is a small value that is almost constant during the period Ts1.
  • period Ta period Te, period Td and period Ts2
  • it increases immediately after the initialization trigger becomes 1, then decreases until the next initialization trigger becomes 1, and converges to an almost constant value. There is.
  • the exponential moving root mean square is a small value that is almost constant during the period Ts1.
  • Ta gradually increases with the increase of the motor torque.
  • Te it converges to an almost constant value.
  • Td it gradually increases as the absolute value of the motor 1009 increases.
  • Ts2 it converges to a small value that is almost constant.
  • the exponential movement skewness is almost 0 in the period Ts1. This is because the motor torque of the period Ts1 has a symmetrical probability distribution. At the time Tr1 at which the period Ta starts, the skewness temporarily becomes a positive value and then becomes a negative value. This is because the motor torque increases from 0, so that old data near 0 and new data after the increase are temporarily mixed, and the probability distribution of the motor torque is disturbed.
  • the exponential movement skewness temporarily becomes a positive value immediately after the time Tr2 in the period Te, and then converges to almost 0.
  • the exponential movement skewness temporarily becomes a positive value immediately after the time Tr3, and then converges to a substantially constant value.
  • the exponential movement skewness temporarily becomes a negative value immediately after the time Tr4, and then converges to almost 0.
  • the exponential movement kurtosis is about 3 in the period Ts1.
  • the value of 3 is a characteristic of the exponential kurtosis with respect to the normal distribution, and indicates that the motor torque has a property close to the normal distribution in the period Ts1.
  • the exponential movement kurtosis becomes a negative value immediately after the time Tr1 at which the period Ta starts, and becomes a positive value thereafter. After that, the exponential kurtosis has converged to a substantially constant value by time Tr2. This is because the motor torque increases from 0, so that old data near 0 and new data after the increase are temporarily mixed, and the probability distribution of the motor torque is disturbed.
  • the exponential movement kurtosis becomes a negative value immediately after the time Tr2 at which the period Te starts, and then becomes a positive value. After that, the exponential kurtosis has converged to a substantially constant value by time Tr3.
  • the exponential movement kurtosis becomes a negative value immediately after the time Tr3 at which the period Td starts, and then becomes a positive value. After that, the exponential kurtosis has converged to a substantially constant value by time Tr4.
  • the exponential movement kurtosis becomes a negative value immediately after the time Tr4 at which the period Ts2 starts, and then becomes a positive value. After that, the exponential movement kurtosis is gradually converging to a constant value by the time Tr5.
  • FIG. 15 is an enlarged view of the time-series waveform of the initialization trigger and the degree of convergence shown in FIG. 14 over a period from time Tr2 to time Tr3.
  • the initialization trigger indicates “1”, and the initialization process is executed.
  • the variable L' j at the jth time point which is one of the internal variables, becomes 0, and the degree of convergence becomes 0.
  • the j-th time point indicates the time point when the initialization trigger indicates "1". The degree of convergence increases after initialization and gradually approaches 1.
  • the state diagnosis unit 4006 diagnoses the state of the mechanical device 1008 based on the feature amount at the time when the convergence degree calculated by the convergence degree calculation unit 4007 using the mathematical formula (85) satisfies a certain condition.
  • the threshold value Cth is set for the degree of convergence.
  • the state diagnosis unit 4006 diagnoses the state of the mechanical device 1008 using the feature amount at the time when the degree of convergence exceeds the threshold value Cth.
  • the threshold value Cth should be set to a value close to 1 in which the degree of convergence is asymptotic.
  • An example of the threshold value Cth is 0.99.
  • Tx the time when the threshold value Cth is exceeded.
  • Tn the period from the time Tr2 to the time Tx in which the degree of convergence does not exceed the threshold value Cth
  • Tc the period from the time Tx when the degree of convergence exceeds the threshold value Cth to Tr3
  • the state diagnosis unit 4006 does not use the feature amount calculated in the period Tn for the diagnosis of the mechanical device 1008, but uses the feature amount calculated in the period Tc for the diagnosis of the mechanical device 1008. According to such a configuration, it is possible to use the feature amount when the calculation result is sufficiently converged when the internal variable is initialized. Thereby, in addition to the effect of the third embodiment, the state of the mechanical device can be estimated with higher accuracy.
  • the information processing apparatus further includes a convergence degree calculation unit with respect to the configuration of the third embodiment.
  • the convergence degree calculation unit calculates the degree of convergence, which is an index that quantitatively indicates the degree of convergence in the calculation of the feature quantity, based on the internal variables.
  • FIG. 16 is a diagram showing a configuration example of the information processing system 100D including the information processing apparatus 5000 according to the fifth embodiment.
  • the state diagnosis unit 1006 is replaced with the state diagnosis unit 5006.
  • Other configurations are the same as or equivalent to the information processing system 100 shown in FIG.
  • the same or equivalent components are designated by the same reference numerals, and duplicate description will be omitted.
  • an example of applying the state diagnosis unit 5006 to the information processing apparatus 1000 shown in FIG. 1 will be described, but the information processing apparatus 2000, 3000, 4000 shown in FIGS. 9, 11 or 13 will be described. It may be applied to any of them.
  • FIG. 17 is a diagram showing a configuration example of the state diagnosis unit 5006 in the fifth embodiment.
  • the state diagnosis unit 5006 includes a state quantity acquisition unit 5401, a learning unit 5402, an abnormality degree calculation unit 5403, and a decision-making unit 5404.
  • the state quantity acquisition unit 5401 acquires a state quantity including a feature quantity.
  • the learning unit 5402 learns the relationship between the state of the mechanical device 1008 and the feature amount based on the state amount of the mechanical device 1008 in the normal state.
  • the abnormality degree calculation unit 5403 calculates the abnormality degree, which is an index quantitatively indicating the abnormality degree of the mechanical device 1008, based on the learning result learned by the learning unit 5402.
  • the decision-making unit 5404 determines the diagnosis result of diagnosing the state of the mechanical device 1008 based on the abnormality degree calculated by the abnormality degree calculation unit 5403.
  • the state quantity acquisition unit 5401 acquires information including a set value given to the mechanical device 1008 as a state quantity together with the feature quantity.
  • a set value given to the mechanical device 1008 As an example of the feature amount, it is conceivable to use the skewness calculated sequentially by the above method. Further, it is conceivable to use the set value of the motor speed as the information of the mechanical device 1008 used as one of the state quantities.
  • the device control unit 1099 is set so that the maximum motor speed is 500 [r / min]. In the present embodiment, it is assumed that the set value of the motor speed is 500 [r / min].
  • the input information to the state quantity acquisition unit 5401 is the set value of the skewness and the motor speed
  • the input information to the state quantity acquisition unit 5401 may be 3 or more.
  • FIG. 18 is a diagram showing the relationship between the feature amount and the set value of the motor speed in the fifth embodiment.
  • a plurality of sets of the skewness which is an example of the feature amount, and the set value of the motor speed are plotted as one data.
  • the horizontal axis is the set value of the motor speed, and the vertical axis is the skewness.
  • the feature quantity may be plural, and the state quantity may be displayed by a three-dimensional concept. Further, a state quantity other than the motor speed may be used.
  • the black circles are normal data in which the mechanical device 1008 is in a normal state.
  • the white circles are abnormality data when some abnormality occurs in the mechanical device 1008. As shown in FIG.
  • the mechanical device 1008 obtains a large number of state quantities in a normal state as normal data. Looking at FIG. 18, it can be seen that the normal data tends to have a larger feature amount as the set value of the motor speed is larger. On the other hand, the abnormal data has a relatively large feature amount even though the set value of the motor speed is relatively small.
  • a method of calculating the abnormality degree using principal component analysis is known.
  • Principal component analysis is a method of re-axising the state quantity, which is multidimensional data, in order from the direction of the largest variance.
  • the learning unit 5402 using the principal component analysis calculates the eigenvalues and eigenvectors of the variance-covariance matrix of the normal data obtained in advance, and linearly maps the original state quantity to the space of the principal component. Assuming that the original state quantity is the vector x and the vector after mapping is y, the mapping of the vector x can be expressed by the following mathematical formula (86).
  • the learning unit 5402 determines the expression matrix A, which is a matrix representing this linear map, based on the normal data.
  • the learning unit 5402 obtains a plurality of eigenvalues and eigenvectors related to the variance-covariance matrix obtained from a plurality of prepared normal data. With respect to the obtained plurality of eigenvalues, a plurality of corresponding eigenvectors are arranged in descending order of the eigenvalues. Then, the matrix formed by the plurality of eigenvectors created side by side becomes the expression matrix A.
  • an eigenvector of 2 rows and 1 column is created using two state quantities.
  • two eigenvectors of 2 rows and 1 column are arranged side by side to create a 2 row and 2 column expression matrix.
  • the dimension of the vector after mapping is made smaller than the dimension of the vector of the original state quantity. May be good.
  • FIG. 19 is a diagram showing the results of principal component analysis in the fifth embodiment.
  • the data of FIG. 18 is reprinted, and the axis of the first principal component having the largest eigenvalue and the axis of the second principal component having the second largest eigenvalue are shown.
  • Eigenvalues are obtained as a result of principal component analysis.
  • FIG. 19 in order to make the result of the principal component analysis easy to understand, according to the eigenvalues corresponding to the variance of each principal component, centering on the position where the axis of the first principal component and the axis of the second principal component intersect. , An elliptical confidence interval is shown.
  • the confidence interval is the interval that is expected to fit in the distribution of the data with a certain probability. For example, there are 90% confidence intervals that are expected to fit 90% of normal data, 99% confidence intervals that are expected to fit 99% of normal data, and so on. Confidence intervals can be calculated using common statistical methods.
  • the normal data is located inside the 99% confidence interval, while the abnormal data is located outside the 99% confidence interval. If the mechanical device 1008 is normal, this abnormal data is judged to be data that is considered to occur only with a probability of about 1%. Therefore, it is determined that there is a possibility that an abnormality has occurred in the mechanical device 1008. On the other hand, if a threshold value is simply set for the feature amount without using principal component analysis, the abnormal data may fall within the distribution of normal data, making diagnosis difficult.
  • the abnormality degree calculation unit 5403 of the fifth embodiment calculates the abnormality degree which is an index which quantitatively indicates the degree of abnormality of the mechanical device 1008.
  • An index called the T2 statistic is used as an example of the degree of anomaly.
  • the T2 statistic is an index in which the data after mapping by principal component analysis standardizes the distance from the center coordinates of the axes of each principal component by the variance for each axis. The distance from the center coordinates of the axis of each principal component is also called "Mahalanobis distance".
  • FIG. 20 is a diagram showing the degree of abnormality of the principal component analysis in the fifth embodiment.
  • FIG. 20 shows a diagram when the normal data, the abnormal data, the 90% confidence interval, and the 99% confidence interval shown in FIG. 19 are standardized using the above-mentioned eigenvalues. Further, in FIG. 20, the distance between the point where the axis of the first principal component and the axis of the second principal component intersect and the abnormality data is shown as the degree of abnormality. Further, the 90% confidence interval and the 99% confidence interval, which were elliptical in FIG. 19, are circular in FIG. 20. As shown in FIG. 20, if the Mahalanobis distance is illustrated using the concept, it is possible to determine the degree of abnormality in consideration of the variation in the data at the normal time.
  • the abnormality degree calculation unit 5403 described a method of calculating the abnormality degree by using the T2 statistic after mapping to the principal component obtained by the principal component analysis, but the present invention is not limited to this method.
  • an index called the Q statistic may be used.
  • unsupervised learning as a method that can learn normal data, remember its distribution, and calculate the degree of abnormality for new data in the same way as principal component analysis.
  • a one-class support vector machine One class Support Vector Machine
  • a Mahalanobis Taguchi Method a self-organizing map (Self-Organizing Maps), or the like
  • the degree of abnormality may be calculated after performing preprocessing such as normalization processing on the input state quantity.
  • the decision-making unit 5404 of the fifth embodiment determines the result of diagnosing the state of the mechanical device 1008 using the degree of abnormality calculated by the above method.
  • a procedure for the decision-making unit 5404 to diagnose the state of the mechanical device 1008 using the feature amount will be described with reference to FIG.
  • FIG. 21 is a diagram showing an example of the secular variation of the degree of abnormality in the fifth embodiment.
  • FIG. 21 shows the result of periodically plotting the degree of anomaly calculated by the above method. Regular means, for example, monthly.
  • the horizontal axis of FIG. 21 is time, which represents the operating time of the mechanical device 1008.
  • the decision-making unit 5404 provides a threshold value Fthb for the degree of abnormality, and when the degree of abnormality exceeds the threshold value Fthb, it diagnoses that an abnormality has occurred in the mechanical device 1008.
  • the time Ttb0 is the time when the operation is started.
  • the degree of abnormality from the time Ttb0 to the time before the time Ttb1 keeps a value smaller than the threshold value Fthb, although there is some variation.
  • the decision-making unit 5404 makes a diagnosis that the mechanical device 1008 is normal, and outputs the diagnosis result.
  • the state diagnosis unit 5006 makes a diagnosis that an abnormality has occurred in the mechanical device 1008 at the time Ttb1, and outputs the diagnosis result.
  • the threshold value Fthb there are various ways to determine the threshold value Fthb. For example, it is conceivable to set the threshold value Fthb as the degree of abnormality corresponding to the 99% confidence interval of normal data. Further, the degree of abnormality when an abnormality occurs in another mechanical device in the past may be set as the threshold value Fthb. Further, the degree of abnormality may be determined with reference to the degree of abnormality immediately after the operation of the mechanical device 1008 is started. Further, the threshold value Fthb may be set dynamically. For example, while operating a plurality of mechanical devices, the degree of abnormality may be calculated individually, and the threshold value Fthb may be determined periodically in consideration of the variation in the degree of abnormality for each device. When the information processing device is configured in this way, when one of the plurality of mechanical devices fails, the degree of abnormality of the failed device becomes larger than the degree of abnormality of the other devices, so that it is possible to detect the abnormality. Become.
  • the method is not limited to this method.
  • a statistical method such as performing a test based on the distribution of the degree of abnormality for a certain period may be adopted.
  • the information processing apparatus diagnoses the state of the mechanical device based on the degree of abnormality, which is an index quantitatively indicating the degree of abnormality of the mechanical device. As a result, it is possible to obtain the effect of the other embodiments and to estimate the state of the mechanical device with high accuracy.
  • FIG. 22 is a diagram showing a configuration example of the information processing system 100E including the information processing apparatus 6000 according to the sixth embodiment.
  • the state diagnosis unit 1006 is replaced with the state diagnosis unit 6006 in the configuration of the information processing system 100 shown in FIG.
  • Other configurations are the same as or equivalent to the information processing system 100 shown in FIG.
  • the same or equivalent components are designated by the same reference numerals, and duplicate description will be omitted.
  • an example of applying the state diagnosis unit 6006 to the information processing apparatus 1000 shown in FIG. 1 will be described, but the information processing apparatus 2000, 3000, 4000 shown in FIGS. 9, 11 or 13 will be described. It may be applied to any of them.
  • FIG. 23 is a diagram showing a configuration example of the state diagnosis unit 6006 in the sixth embodiment.
  • the state diagnosis unit 6006 includes a state quantity acquisition unit 6401, a learning unit 6402, a state estimation unit 6403, and a decision-making unit 6404.
  • the state quantity acquisition unit 6401 acquires a state quantity including a feature quantity.
  • the learning unit 6402 learns the relationship between the state of the mechanical device 1008 and the state amount based on the learning data that associates the state of the mechanical device 1008 with the state amount acquired by the state amount acquisition unit 6401.
  • the state estimation unit 6403 estimates the state of the mechanical device 1008 based on the state amount input to the learning unit 6402 and the learning result by the learning unit 6402.
  • the decision-making unit 6404 determines the diagnostic result of diagnosing the state of the mechanical device 1008 based on the result estimated by the state estimation unit 6403.
  • the state quantity acquisition unit 6401 acquires information including a set value given to the mechanical device 1008 as a state quantity together with the feature quantity.
  • a set value given to the mechanical device 1008 As an example of the feature amount, it is conceivable to use the skewness and kurtosis calculated sequentially by the above method. Further, it is conceivable to use the above-mentioned set value of the motor speed as the information of the mechanical device 1008 used as one of the state quantities.
  • three cases where the input information to the state quantity acquisition unit 6401 is the set value of skewness, kurtosis, and motor speed are exemplified, but the present invention is not limited to this example.
  • the input information to the state quantity acquisition unit 6401 may be 4 or more.
  • the learning unit 6402 and the state estimation unit 6403 may learn the relationship between the state quantity and the state of the mechanical device 1008 by so-called supervised learning, for example, according to a neural network model.
  • supervised learning by giving a large amount of learning data, which is a set of data of a certain input (state amount) and result (label), to the learning unit, a model that learns the features of those data and estimates the result from the input. Is called supervised learning.
  • a neural network is composed of an input layer composed of a plurality of neurons, an intermediate layer (hidden layer) composed of a plurality of neurons, and an output layer composed of a plurality of neurons.
  • the intermediate layer may be one layer or two or more layers.
  • FIG. 24 is a diagram showing an example of the structure of the state estimation unit 6403 in the sixth embodiment.
  • the input is a state quantity.
  • the neural network of FIG. 24 has 3 inputs and 3 layers.
  • the value obtained by multiplying the input value by the weight W1 composed of W11 to W16 is input to the intermediate layer composed of Y1 and Y2. Will be done.
  • the value obtained by multiplying the input value of the intermediate layer by the weight W2 composed of W21 to W26 is output from the output layer composed of Z1 to Z3. This output result changes depending on the value of the weight W1 and the value of the weight W2.
  • the state quantities 1 to 3 are input to the input layers X1 to X3, and the probability of each state output from the output layers Z1 to Z3 matches the label of the training data.
  • the weights W1 and W2 are adjusted so as to do so.
  • the information processing apparatus 6000 equipped with the learned state estimation unit 6403 that has executed the learning process in the sixth embodiment may be configured.
  • the trained state estimation unit 6403 may be composed of trained data, trained programs, or a combination thereof. By using the learned state estimation unit 6403, the result of learning using another information processing device can be used, so that an information processing device 6000 capable of realizing diagnosis without performing new learning is provided. be able to.
  • FIG. 24 shows an example in which the skewness is input to the input layer X1, the kurtosis is input to the input layer X2, and the set value of the motor speed is input to the input layer X3.
  • the strain degree is an example of the state quantity 1 input to the input layer X1
  • the sharpness is an example of the state quantity 2 input to the input layer X2
  • the set value of the motor speed is input to the input layer X3.
  • the probability that the ball screw shaft 1224 of the mechanical device 1008 has failed is output as the probability of the state 1.
  • the probability that the coupling 1220 of the mechanical device 1008 has failed is output as the probability of the state 2.
  • the probability that the mechanical device 1008 is in a normal state is output as the probability of the state 3.
  • a softmax function softmax function
  • this function makes the estimation result easier to understand.
  • the neural network shown in FIG. 24 learns the relationship between the state of the mechanical device 1008 and the state quantities 1 to 3 by supervised learning according to a data set created based on the learning data input to the learning unit 6402. .
  • the information indicating the true state of the mechanical device 1008 called the label is correctly obtained, and the output of any one of the states 1 to 3 is 1 (probability is 100%), and the other outputs are. It is better to learn so that it becomes 0 (probability is 0%).
  • the configuration of the learning unit 6402 and the state estimation unit 6403 does not have to be limited to the neural network.
  • there are many known methods such as k-nearest neighbor method, binary tree search, support vector machine, linear regression, logistic regression, etc., in which the relationship between the state quantity and the label is learned and the state is estimated from the state quantity. Therefore, the learning unit 6402 and the state estimation unit 6403 may be configured by applying this method.
  • the decision-making unit 6404 determines the diagnosis result of diagnosing the state of the mechanical device 1008 based on the estimation result of the state estimation unit 6403. For example, when the mechanical device 1008 is most likely to be normal, it is preferable to output a diagnostic result indicating that the mechanical device 1008 is normal as a diagnostic result. Further, when the probability that the ball screw shaft 1224 has failed or the probability that the coupling 1220 has failed is the highest, it is preferable to output the diagnosis result indicating the failure location. Further, when the sum of the probability that the ball screw shaft 1224 has failed and the probability that the coupling 1220 has failed is larger than the normal probability, a result indicating that any part has failed is output. It is better to do it.
  • the information processing apparatus learns the relationship between the state of the mechanical device and the state quantity based on the learning data that links the relationship between the state of the mechanical device and the state quantity.
  • the state of the mechanical device is estimated based on the learning result and the state quantity used at the time of learning.
  • the configuration shown in the above embodiments is an example, and can be combined with another known technique, can be combined with each other, and does not deviate from the gist. It is also possible to omit or change a part of the configuration.
  • 100,100A, 100B, 100C, 100D, 100E Information processing system 1000,2000,3000,4000,5000,6000 Information processing device, 1001,2001,3001 Sensor data acquisition unit, 1002 Internal variable holding unit, 1003,2003 4003 Internal variable calculation unit, 1004, 2004, 4004 Feature quantity calculation unit, 1005, 3005 Initialization processing unit, 1006,400,500,6006 State diagnosis unit, 1008 mechanical device, 1009 motor, 1010 sensor, 1099, 2099 device control Part, 1210 ball screw, 1212 movable part, 1213 guide, 1220 coupling, 1224 ball screw shaft, 1230 servo motor, 1231 servo motor shaft, 1232 current sensor, 1233 encoder, 1240 driver, 1250, 1280 indicator, 1260 PLC, 1270 PC, 1291 processor, 1292 memory, 1293 processing circuit, 4007 convergence degree calculation unit, 5401, 6401 state quantity acquisition unit, 5402, 6402 learning unit, 5403 abnormality degree calculation unit, 6403 state estimation unit, 5404, 6404 decision unit ..

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)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
  • Testing And Monitoring For Control Systems (AREA)
  • Hardware Redundancy (AREA)

Abstract

情報処理装置（１０００）において、内部変数保持部（１００２）は、センサデータ取得部（１００１）によって取得されたセンサデータに基づいて時系列に順次計算される、Ｎよりも少ない個数の内部変数を保持する。内部変数計算部（１００３）は、第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する内部変数を、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算する。特徴量計算部（１００４）は、第Ｎ時点の内部変数に基づいて第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量を計算する。状態診断部（１００６）は、特徴量に基づいて機械装置（１００８）の状態を診断する。

Description

情報処理装置及び情報処理方法
　本開示は、機械装置の状態を診断するための情報処理を実施する情報処理装置及び情報処理方法に関する。
　例えば、ボールねじ、減速機などの機械部品が含まれる機械装置では、機械部品が経年劣化し、摩擦の増加、振動の発生、筐体の破損といった様々な異常が生じる。このため、この種の異常を早期に発見して把握する情報処理装置は、効率的な工場での運用に重要とされている。
　機械装置の状態を診断する情報処理装置の例が、下記特許文献１に記載されている。特許文献１には、時系列に取得されるセンサデータを予め設定された期間に渡って保持し、保持したセンサデータから平均値、分散等の特徴量を計算する方法が記載されている。また、特許文献１では、計算した特徴量を予め設定された期間に渡って保持し、保持した特徴量に基づいて歪度等の指標を計算し、計算した値を過去の値と比較することで機械装置の異常を検知することが記載されている。
特開２０２０－３５３７２号公報
　しかしながら、特許文献１に記載の情報処理装置においては、センサの信号によって得られるセンサデータを設定期間に渡って保持する必要があり、膨大な記憶容量が要求されるという課題がある。例えば、センサデータを得るためのサンプリング周波数が２ｋＨｚであれば、１秒間に２千点の値が得られる。ここで、センサデータを計算機で扱うためのデータ型をｉｎｔ型４バイトとすると、例えば６０秒間のセンサデータを得るためには４８０キロバイト（＝４バイト×２ｋＨｚ×６０秒間）のメモリが必要となる。産業用に使われる多くのマイクロコントローラのＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）が数キロから数メガバイト程度であることから、情報を処理可能な装置が限定されてしまうという問題があった。
　また、センサデータをマイクロコントローラではなくサーバ装置のストレージに保持して別途オフラインで処理する方法も考えられる。しかしながら、このようなストレージは、サーバ装置を設置するスペース、サーバ装置を管理するための管理コストなどの観点から、工場に数台程度しか用意できないのが現状である。このため、数台のストレージでは、日々蓄積されていく膨大なデータを処理しきれないという問題がある。例えば、工場に前述のセンサが１０個備えられた機械装置が１０台あり、それらが２４時間稼働する場合を考える。この場合、一ヶ月毎に２テラバイト以上（≒１０個×１０台×４バイト×２［ｋＨｚ］×３６００秒×２４時間×３０日＝２．０７３６テラバイト）のデータが追加されることになり、サーバ装置の計算機にとって大きな負担となる。
　本開示は、上記に鑑みてなされたものであって、計算能力の高い計算機及び記憶容量の大きな計算機を用いることなく、機械装置の状態を診断するための情報処理を実施できる情報処理装置を得ることを目的とする。
　上述した課題を解決し、目的を達成するため、本開示に係る情報処理装置は、センサデータ取得部と、内部変数保持部と、内部変数計算部と、特徴量計算部と、状態診断部とを備える。センサデータ取得部は、センサによって計測された機械装置の物理量の計測値を取得し、当該計測値のうちの第１時点から第Ｎ時点（Ｎは２以上の整数）までの計測値をセンサデータとして保持する。内部変数保持部は、センサデータに基づいて時系列に順次計算される、Ｎよりも少ない個数の内部変数を保持する。内部変数計算部は、第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する内部変数を、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算する。特徴量計算部は、第Ｎ時点の内部変数に基づいて第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量を計算する。状態診断部は、特徴量に基づいて機械装置の状態を診断する。
　本開示に係る情報処理装置によれば、計算能力の高い計算機及び記憶容量の大きな計算機を用いることなく、機械装置の状態を診断するための情報処理を実施できるという効果を奏する。
実施の形態１における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態１における機械装置及びその周辺装置のハードウェア構成例を示す図 図２に示すドライバが備える処理回路をプロセッサ及びメモリで構成する場合の構成例を示す図 図２に示すドライバが備える処理回路を専用のハードウェアで構成する場合の構成例を示す図 実施の形態１におけるモータ速度及びモータトルクの時系列の波形を示す図 実施の形態１におけるモータトルク及び各種の特徴量の時系列の波形を示す図 実施の形態１における特徴量の経年変化の様子の一例を示す図 実施の形態１における情報処理方法の説明に供するフローチャート 実施の形態２における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態２におけるモータトルク及び各種の特徴量の時系列の波形を示す図 実施の形態３における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態３におけるモータトルク及び各種の特徴量の時系列の波形を示す図 実施の形態４における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態４におけるモータトルク及び各種の特徴量の時系列の波形を示す図 図１４に示す初期化トリガ及び収束度の時系列の波形を時刻Ｔｒ２から時刻Ｔｒ３までの期間に渡って拡大した図 実施の形態５における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態５における状態診断部の構成例を示す図 実施の形態５における特徴量とモータ速度の設定値との関係を示す図 実施の形態５における主成分分析の結果を示す図 実施の形態５における主成分分析の異常度を示す図 実施の形態５における異常度の経年変化の様子の一例を示す図 実施の形態６における情報処理装置を含む情報処理システムの構成例を示す図 実施の形態６における状態診断部の構成例を示す図 実施の形態６における状態推定部の構造の例を示す図
　以下に、本開示の実施の形態における情報処理装置及び情報処理方法を図面に基づいて詳細に説明する。
実施の形態１．
　図１は、実施の形態１における情報処理装置１０００を含む情報処理システム１００の構成例を示す図である。情報処理システム１００は、図１に示すように、情報処理装置１０００と、機械装置１００８と、センサ１０１０と、モータ１００９と、装置制御部１０９９とを備える。モータ１００９は、機械装置１００８に駆動力を付与することで機械装置１００８を駆動する。機械装置１００８がモータ１００９によって駆動される際、センサ１０１０は機械装置１００８の物理量を計測する。機械装置１００８が計測する物理量の例は、位置、速度、加速度、動作指令、電流、電圧、トルク、力、圧力、音声、又は光量である。なお、本実施の形態では、モータ速度及びモータトルクを一例として説明する。モータ速度はモータ１００９の回転速度であり、モータトルクはモータ１００９が発生するトルクである。
　センサ１０１０は、物理量の計測値を含むセンサ信号を情報処理装置１０００及び装置制御部１０９９に出力する。装置制御部１０９９は、センサ信号に基づいて、モータ１００９を制御するための制御信号を決定する。モータ１００９は、装置制御部１０９９から出力される制御信号によって制御される。
　情報処理装置１０００は、センサデータ取得部１００１と、内部変数保持部１００２と、内部変数計算部１００３と、特徴量計算部１００４と、初期化処理部１００５と、状態診断部１００６とを備える。
　センサデータ取得部１００１は、センサ１０１０によって計測された機械装置１００８の物理量の計測値を取得する。前述したように、物理量の計測値は、センサ１０１０が送信するセンサ信号に含まれている。また、センサデータ取得部１００１は、取得した計測値のうちの第１時点から第Ｎ時点までの計測値をセンサデータとして保持する。ここで、Ｎは２以上の整数である。内部変数保持部１００２は、当該Ｎよりも少ない個数の内部変数を一時的に保持する。内部変数は、センサデータに基づいて時系列に順次計算される変数である。内部変数は、特徴量の計算に用いられる。内部変数及び特徴量の詳細は、後述する。
　内部変数計算部１００３は、センサデータ取得部１００１から送信されるセンサデータと、内部変数保持部１００２から送信される内部変数とを受信する。内部変数計算部１００３は、センサデータと内部変数とに基づいて内部変数を更新する。より一般化して説明すると、内部変数計算部１００３は、第ｊ＋１時点に対応する内部変数を、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算することで、内部変数を逐次的に更新していく。ここで、ｊは１からＮ－１までの整数である。即ち、内部変数計算部１００３は、ある時点でのセンサデータと、ある時点よりも１時点前の内部変数に基づいて、当該ある時点における内部変数を計算する。計算された内部変数は、内部変数保持部１００２に送信される。
　特徴量計算部１００４は、センサデータ取得部１００１から送信されるセンサデータと、内部変数保持部１００２から送信される内部変数とを受信する。特徴量計算部１００４は、センサデータ及び内部変数に基づいて特徴量を計算する。より一般化して説明すると、特徴量計算部１００４は、内部変数に基づいて、第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量を計算する。計算された特徴量は、状態診断部１００６に送信される。
　初期化処理部１００５は、初期化処理を実行する。初期化処理は、内部変数保持部１００２が保持する内部変数を初期値に設定する処理である。より一般化して説明すると、初期化処理部１００５は、第１時点の内部変数を、予め設定した最大値と予め設定した最小値との間の値に決定する処理を行う。
　状態診断部１００６は、特徴量に基づいて機械装置１００８の状態を診断する診断処理を行い、診断処理の結果である診断結果を出力する。
　図２は、実施の形態１における機械装置１００８及びその周辺装置のハードウェア構成例を示す図である。図２には、機械装置１００８の構成例として、サーボモータ１２３０を駆動源として用いた機械装置１００８が示されている。サーボモータ１２３０が発生した駆動トルクは、サーボモータ軸１２３１から出力され、カップリング１２２０を介してボールねじ軸１２２４に入力される。ボールねじ１２１０は、回転動作をねじ機構で直動の動作に変換し、可動部１２１２を動作させる。なお、図示は省略するが、可動部１２１２は、異なる機械部品に接続され、可動された機械部品が機械装置１００８の目的に応じて利用される。
　可動部１２１２は、ガイド１２１３で所望の方向へ移動を制限される。ガイド１２１３は、機械装置１００８が精度良く動作できるように可動部１２１２を補助している。サーボモータ１２３０には、図示のように、サーボモータ軸１２３１を予め定めた、位置、速度又はトルクに追従して駆動できるように、回転角度を計測するエンコーダ１２３３と、電流を計測する電流センサ１２３２とが併設されているのが一般的である。ドライバ１２４０は、電流センサ１２３２から得た情報を元にフィードバック制御を行い、駆動に必要な電力をサーボモータ１２３０に供給する。フィードバック制御に必要な計算は、図１の装置制御部１０９９が実施する。なお、本実施の形態では、状態を診断するために用いるセンサデータの例として、モータトルクを例示するが、これに限定されない。機械装置１００８の状態が情報として含まれるものならよく、前述した、位置、速度、加速度、電流、電圧、トルク、力、圧力、音声、光量などの物理量でもよい。また、これらの物理量に代えて、画像情報などでもよい。
　また、図２では、センサ１０１０の例として、電流センサ１２３２及びエンコーダ１２３３を例示しているが、これらに限定されない。他の例として、レーザ変位計、ジャイロセンサ、振動計、加速度センサ、電圧計、トルクセンサ、圧力センサ、マイク、光センサ、カメラなどが例示できる。また、センサ１０１０の取り付け位置は、必ずしもサーボモータ１２３０に近接している必要はなく、機械装置１００８の状態を診断するのに適している位置であればよい。例えば、加速度センサをガイド１２１３の外面等に設置し、センサデータとして加速度を計測してもよい。
　ＰＬＣ（Ｐｒｏｇｒａｍｍａｂｌｅ　Ｌｏｇｉｃ　Ｃｏｎｔｒｏｌｌｅｒ）１２６０は、サーボモータ１２３０の動作命令をドライバ１２４０へ送る。また、必要に応じてＰＣ（Ｐｅｒｓｏｎａｌ　Ｃｏｍｐｕｔｅｒ）１２７０が、用意される場合がある。この場合、ＰＣ１２７０は、ＰＬＣ１２６０へ命令を発信するために使用される。ＰＣ１２７０として、産業用ＰＣ（Ｆａｃｔｏｒｙ　Ａｕｔｏｍａｔｉｏｎ　ＰＣ又はＩｎｄｕｓｔｒｉａｌ　ＰＣ）が用いられる場合もある。
　また、図示のように、ＰＬＣ１２６０の状態をモニタするＰＬＣ用の表示器１２５０と、ＰＣ１２７０の状態をモニタするＰＣ用の表示器１２８０とが用意される場合がある。なお、図示は省略するが、サーボモータ１２３０等の駆動源は、一つの機械装置１００８に複数設けられることが一般的である。このため、ドライバ１２４０も必要に応じて複数用意されることがある。また、単一のＰＬＣ１２６０が統括して、或いは複数のＰＬＣ１２６０が協調して機械装置１００８を動作させる場合がある。これらの構成の場合にも、本実施の形態で示す情報処理装置１０００は、同様に実施可能である。
　図３は、図２に示すドライバ１２４０が備える処理回路をプロセッサ１２９１及びメモリ１２９２で構成する場合の構成例を示す図である。処理回路がプロセッサ１２９１及びメモリ１２９２で構成される場合、ドライバ１２４０の処理回路の各機能は、ソフトウェア、ファームウェア、又はソフトウェアとファームウェアとの組み合わせによって実現される。ソフトウェア又はファームウェアはプログラムとして記述され、メモリ１２９２に格納される。処理回路では、メモリ１２９２に記憶されたプログラムをプロセッサ１２９１が読み出して実行することによって、各機能を実現する。即ち、処理回路は、ドライバ１２４０の処理が結果的に実行されることになるプログラムを格納するためのメモリ１２９２を備える。また、これらのプログラムは、ドライバ１２４０の手順及び方法をコンピュータに実行させるものであるとも言える。
　ここで、プロセッサ１２９１は、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、又はＤＳＰ（Ｄｉｇｉｔａｌ　Ｓｉｇｎａｌ　Ｐｒｏｃｅｓｓｏｒ）と称される演算手段であってもよい。メモリ１２９２は、例えば、ＲＡＭ、ＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ　Ｐｒｏｇｒａｍｍａｂｌｅ　ＲＯＭ）、ＥＥＰＲＯＭ（登録商標）（Ｅｌｅｃｔｒｉｃａｌｌｙ　ＥＰＲＯＭ）などの、不揮発性又は揮発性の半導体メモリとしてもよい。また、メモリ１２９２を、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、又はＤＶＤ（Ｄｉｇｉｔａｌ　Ｖｅｒｓａｔｉｌｅ　Ｄｉｓｃ）などの、記憶手段としてもよい。
　図４は、図２に示すドライバ１２４０が備える処理回路を専用のハードウェアで構成する場合の構成例を示す図である。処理回路が専用のハードウェアで構成される場合、図４に示す処理回路１２９３は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、又はこれらを組み合わせたものとしてもよい。ドライバ１２４０の機能を、機能毎に処理回路１２９３によって実現してもよく、複数の機能をまとめて処理回路１２９３によって実現してもよい。なお、ドライバ１２４０とＰＬＣ１２６０とは、ネットワークを介して接続してもよい。また、ＰＣ１２７０は、クラウドサーバ上に存在してもよい。
　ハードウェア構成の一例は前述の通りであるが、ドライバ１２４０、ＰＬＣ１２６０、ＰＣ１２７０は必須ではなく、本開示に係る情報処理装置を実施するためのデバイスを別途用意し、当該デバイスの内部で実施してもよい。例えば、バッテリ、マイコン、センサ、表示器、通信機能を備えた単独のデバイスを使用し、機械装置１００８の音声をマイクで取得したセンサデータに基づいて、機械装置１００８の状態を推定する構成としてもよい。
　また、表示器類は必須ではなく、ＰＬＣ用の表示器１２５０又はＰＣ用の表示器１２８０で結果を表示する代わりに、ドライバ１２４０、ＰＬＣ１２６０に備えられた既存のＬＥＤなどを使用して、結果を表示してもよい。また、表示器で結果を示さずに、異常が生じたことを診断した際にサーボモータ１２３０の駆動を停止するような構成にしてもよい。
　また、本実施の形態では、モータ速度をサーボモータ１２３０に備えられたエンコーダ１２３３から得ることとして説明したが、これに限定されない。例えば、図１の構成において、モータ１００９へ駆動の命令を下す装置制御部１０９９からの制御信号を用いてモータ速度を得るようにしてもよい。
　また、本実施の形態において、サーボモータ１２３０は回転型のサーボモータとして説明としたが、リニア型のサーボモータ、誘導機モータ、ステッピングモータ、ブラシモータ、超音波モータ等、その他のモータ又は駆動源を用いて実施してもよい。また、ボールねじ１２１０及びカップリング１２２０は、構成品の例であってこれらに限定されない。減速機、ガイド、ベルト、スクリュー、ポンプ、ベアリング、筐体等、その他多様な部品で構成された機械装置に対しても、本開示に係る情報処理装置への適用は可能である。
　次に、実施の形態１における情報処理装置１０００の動作について、図面及び数式を用いて詳細に説明する。まず、サーボモータ１２３０を用いて機械装置１００８を駆動した際に電流センサ１２３２から得られるセンサデータに関して、図５を用いて説明する。図５は、実施の形態１におけるモータ速度及びモータトルクの時系列の波形を示す図である。
　電流センサ１２３２から直接得られる信号は、モータ１００９内を流れる三相の電流を計測した信号である。モータトルクは、三相の電流に適切な変換を施すことで計算することができる。なお、電流センサ１２３２によって検出された三相の電流値を、センサデータとして診断に用いてもよい。また、エンコーダ１２３３から得られる信号は、モータ１００９の回転角度を表す位置情報である。従って、位置情報に対して数値微分等の処理を施せば、モータ１００９の回転速度であるモータ速度が得られる。このため、電流センサ１２３２から得られる信号に代えて、エンコーダ１２３３から得られる信号をセンサデータとして診断に用いてもよい。
　図５（ａ）、（ｂ）の横軸は、時間を表している。図５（ａ）にはモータ速度の時系列の波形が示され、図５（ｂ）にはモータトルクの時系列の波形が示されている。これらは、モータ１００９が停止している状態から位置決めと呼ばれる単発の駆動を実施した際の波形である。また、図５（ｂ）の波形の取得周期であるサンプリング周期は０．５ミリ秒であり、取得時間は６秒である。従って、センサデータの点数であるデータ点数Ｎは、１２０００点（＝６０００ミリ秒÷０．５ミリ秒）である。なお、ここで示すサンプリング周期、取得時間及びデータ点数Ｎは一例であり、これらの数値に限定されるものではない。また、本例では、位置決めの回数は１回としているが、位置決めの回数は複数回でもよい。また、ここでは、モータ１００９を位置決めする場合の動作について説明するが、これに限定されない。本実施の形態は、位置決め以外の制御、例えば速度制御又はトルク制御への適用も可能である。
　図５（ａ）のモータ速度に関して説明する。まず、時刻Ｔｒ０から時刻Ｔｒ１の間は、モータ１００９は停止しており、モータ速度は０［ｒ／ｍｉｎ］である。この期間を「Ｔｓ１」と表記する。続いて、時刻Ｔｒ１から時刻Ｔｒ２の間にモータ１００９は加速し、モータ速度は５００［ｒ／ｍｉｎ］まで上昇する。この期間を「Ｔａ」と表記する。続いて、時刻Ｔｒ２から時刻Ｔｒ３の間は、モータ速度は一定で５００［ｒ／ｍｉｎ］のままである。この期間を「Ｔｅ」と表記する。続いて、時刻Ｔｒ３から時刻Ｔｒ４の間にモータ１００９は減速し、モータ速度は０［ｒ／ｍｉｎ］まで下降する。この期間を「Ｔｄ」と表記する。続いて、時刻Ｔｒ４から時刻Ｔｒ５の間は、モータ１００９は停止しており、モータ速度０［ｒ／ｍｉｎ］のままである。この期間を「Ｔｓ２」と表記する。以上の動作を図示すると、図５（ａ）に示されるような台形状の波形となる。なお、このような台形状の波形は一例であり、これに限定されない。モータ１００９の駆動に伴うセンサデータから機械装置１００８の状態が得られる動作であればよく、モータ速度の波形の形状はどのようなものでもよい。
　次に、図５（ｂ）のモータトルクに関して説明する。期間Ｔｓ１において、モータ１００９は停止しており、モータ１００９の動作に必要なモータトルクはほぼ０［Ｎｍ］である。続く期間Ｔａでは、機械装置１００８を加速させるためのトルクが必要となり、期間Ｔｓ１よりも大きなモータトルクが発生する。なお、粘性摩擦等の影響により機械装置１００８に粘性摩擦が存在すると、図５のように、速度に応じて必要なトルクが徐々に大きくなる場合もある。続く期間Ｔｅでは、機械装置１００８の速度は変わらず、ほぼ一定のモータトルクが発生している。続く期間Ｔｄでは、機械装置１００８を減速させるためのトルクが必要となり、期間Ｔｅのモータトルクから見て負方向にモータトルクが発生する。最後の期間Ｔｓ２において、モータ１００９は停止しており、モータ１００９の動作に必要なモータトルクは、ほぼ０［Ｎｍ］である。
　以上に説明したモータトルクは、ほぼ理想的な機械装置１００８の動作に必要なトルクである。但し、実際には、電気的な振動、又は機械的な振動、摩擦等の影響により、図５（ｂ）に示すような、大なり小なりのノイズがモータトルクに発生する。機械装置１００８の部品であるボールねじ軸１２２４、ガイド１２１３、カップリング１２２０等が劣化した場合、前述の摩擦又は振動の大きさが変化し、モータトルクを含むセンサデータに摩擦又は振動の影響が現れることが知られている。そのため、センサデータを統計的な手法で解析し、特徴量と呼ばれる幾つかの指標を計算して装置の異常を検出することが行われる。
　センサデータ取得部１００１は、センサ１０１０から時系列に逐次生成されるセンサ信号を取得し、デジタル化されたセンサデータを生成する。生成されたセンサデータは、内部変数計算部１００３及び特徴量計算部１００４のうちの少なくとも一つに渡される。センサデータ取得部１００１は、センサデータを生成する際、必要に応じて機械装置１００８の状態に関係のないノイズを除去するフィルタ処理を実行してもよい。
　次に、内部変数保持部１００２、内部変数計算部１００３、特徴量計算部１００４及び初期化処理部１００５の機能及び動作について、幾つかの種類の特徴量を例にとり、詳細に説明する。
　まず、時刻及びセンサデータの数式中の表記に関して説明する。特徴量計算に用いるセンサデータの取得を開始した時点を「第１時点」と呼ぶ。センサデータは、第１時点から第Ｎ時点までに、合計でＮ個得られる。ここで、Ｎは２以上の整数である。また、第ｊ番目のセンサデータを取得した時点を「第ｊ時点」と呼ぶ。ここで、ｊは１からＮ－１までの整数である。第ｊ時点で得られたセンサデータを「ｘｊ」と記す。以降、説明を簡単にするため、センサデータは、予め定めた時間間隔毎に取得されることとする。なお、言うまでもないが、センサデータは、不定期に取得されてもよい。
　内部変数保持部１００２は、１時点前の内部変数を保持する。内部変数は１種類とは限らず、特徴量に応じて複数種類保持する場合もある。保持する内部変数の種類が、特徴量の計算の基となるセンサデータの時系列の数よりも小さいほど、メモリの削減効果が高い。内部変数保持部１００２は、特徴量それ自体を内部変数の一つの種類として保持してもよい。
　初期化処理部１００５は、保持した内部変数を初期値に設定する初期化処理を実行する。初期値は、初期化処理の際に設定する値である。初期化処理部１００５は、例えば情報処理装置１０００の電源が投入され、センサデータ取得部１００１が最初にセンサデータを取得した時刻である第１時点で初期化処理を実施する。初期値は、後述するように０とするのが簡単であるが、必ずしも０でなくともよい。例えば、適切な初期値の付近に上限値と下限値を定め、初期値を上限値と下限値との間に適宜決定してもよい。なお、後述する一部の特徴量は、内部変数及び特徴量の初期値を極端に小さい値にすると計算が安定しない場合がある。この場合、初期値を０以外の値にすることで計算が安定化する。また、センサデータに対する定常的な値が０でない特徴量も存在する。そのような特徴量を用いる場合、想定される特徴量の値を初期値とすることで、計算の収束を早めることができる。例えば、本実施の形態で示す特徴量の一つである尖度は、正規分布に近い性質を持つセンサデータに対して３程度の値になることが知られている。そのため、尖度の初期値を３にすることで、多くの種類のセンサに対して収束の早い計算が行える。以降、特徴量の具体的な例として、平均値、分散、標準偏差、二乗平均平方根、歪度、尖度、最大値、最小値、ピーク値、ピークピーク値及び波高率を用いたときの、各特徴量の逐次的な計算方法を、その導出手順と共に示す。なお、以降、各特徴量の逐次的な計算を、単に「逐次計算」と呼ぶ場合がある。
　まず、特徴量の最も単純な一例として、平均値ｍの逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する平均値ｍｊは、以下の公知の数式（１）で定義される。この数式（１）は、算術平均値又は相加平均値とも呼ばれる。
Figure JPOXMLDOC01-appb-M000001
　前述のように、第ｉ時点のセンサデータはｘｊである。従って、第３時点の平均値ｍ３は、以下の数式（２）のように表せる。
Figure JPOXMLDOC01-appb-M000002
　ここで、第１時点から第３時点までの時系列のセンサデータをメモリ１２９２には記憶せず、第３時点の情報を幾つかの内部変数として保持する。そして、第３時点の内部変数と第４時点のセンサデータとに基づいて、第４時点の平均値を求める方法を考える。まず、第４時点の平均値ｍ４は、定義式である数式（１）に従って、以下の数式（３）のように表せる。
Figure JPOXMLDOC01-appb-M000003
　上記の数式（３）を用いて平均値ｍ４を算出する場合、第１から第４時点までの時系列のセンサデータを用いており、過去のサンプル数分のセンサデータを保持するメモリが必要となる。しかし、第３時点での平均値ｍ３と、第３時点までのセンサデータの点数を表す変数Ｌ３＝３を内部変数として保持すれば、第４時点の平均値ｍ４は、以下の数式（４）で表せる。
Figure JPOXMLDOC01-appb-M000004
　上記の数式（４）に、数式（２）及びＬ３＝３を代入すると、平均値ｍ４を定義式で求めた数式（３）の結果と一致する。
　一般化して考えると、第ｊ＋１時点の平均値ｍｊ＋１は、以下の数式（５）のように表せる。
Figure JPOXMLDOC01-appb-M000005
　即ち、第ｊ＋１時点の平均値ｍｊ＋１は、以下の数式（６）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000006
　ここで、第ｊ＋１時点の変数Ｌｊ＋１は、以下の数式（７）で逐次的に求まる。
Figure JPOXMLDOC01-appb-M000007
　上記の数式（７）で第ｊ時点の変数を保持することは、逐次計算中の時刻ｊの情報を保持していることと等しい。センサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第ｊ＋１時点の平均値ｍｊ＋１を求めるために第ｊ時点で保持すべき内部変数は、平均値ｍｊ及び変数Ｌｊである。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ及び平均値ｍの初期値Ｌ１，ｍ１は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値に応じて決定するのがよい。例えば、特徴量に平均値を用いる場合、変数Ｌの初期値Ｌ１＝０、平均値ｍの初期値ｍ１＝ｘ１とすると、結果が定義式である数式（１）と一致する。
　次に、特徴量の一例として、分散ｖの逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する分散ｖｊは、以下の公知の数式（８）で定義される。
Figure JPOXMLDOC01-appb-M000008
　なお、上記の数式（８）の代わりに、分母をｊ－１とする不偏分散を用いることもあるが、その場合も同様の手順で導出可能である。
　ここで、センサデータを確率変数と見なすと、その分散ｖは、以下の数式（９）で表せるということが知られている。
Figure JPOXMLDOC01-appb-M000009
　上記数式（９）において、記号“Ｅ［］”は、角括弧内の確率変数の期待値を表している。従って、上記数式（９）を用いて上記数式（８）に示される第ｊ時点の分散ｖｊは、以下の数式（１０）で表せる。
Figure JPOXMLDOC01-appb-M000010
　従って、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する分散ｖｊ＋１は、以下の数式（１１）で表せる。
Figure JPOXMLDOC01-appb-M000011
　また、上記の数式（１１）を変形すると、分散ｖｊ＋１は、以下の数式（１２）で表せる。
Figure JPOXMLDOC01-appb-M000012
　そして、上記の数式（１０）と数式（１２）とから、以下の数式（１３）が導かれる。
Figure JPOXMLDOC01-appb-M000013
　従って、上記の数式（１３）より、第ｊ＋１時点の分散ｖｊ＋１は、以下の数式（１４）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000014
　上記のように、第ｊ＋１時点の変数Ｌｊ＋１は、前述の数式（７）で逐次的に求められる。また、第ｊ＋１時点の平均値ｍｊ＋１は、前述の数式（６）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第ｊ＋１時点の分散ｖｊ＋１を求めるために第ｊ時点で保持すべき内部変数は、変数Ｌｊ、分散ｖｊ及び平均値ｍｊである。なお、分散ｖｊに代えて標準偏差ｓｊを保持してもよい。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ及び分散ｖの初期値Ｌ１，ｖ１は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。ただし、分散ｖの初期値ｖ１を０に近い値に設定すると、後述の別の特徴量の逐次計算で分散ｖが分母に現れる場合に計算が安定しない。このため、０よりもある程度大きい値を初期値に用いるとよい場合もある。
　次に、特徴量の一例として、標準偏差ｓの逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する標準偏差ｓｊは、以下の公知の数式（１５）で定義される。
Figure JPOXMLDOC01-appb-M000015
　即ち、標準偏差ｓｊは、前述の手順で逐次的に求めた分散ｖｊから直ちに求まる。そのため、標準偏差ｓｊを逐次的に計算するために保持すべき内部変数は、分散ｖｊと同じである。
　次に、特徴量の一例として、二乗平均平方根ｒの逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する二乗平均平方根ｒｊは、以下の公知の数式（１６）で定義される。これは、ＲＭＳ（Ｒｏｏｔ　Ｍｅａｎ　Ｓｑｕａｒｅ：実効値）とも呼ばれる。
Figure JPOXMLDOC01-appb-M000016
　従って、第１時点から第ｊ＋１時点までのセンサデータに対する二乗平均平方根ｒｊ＋１は、以下の数式（１７）のように表せる。
Figure JPOXMLDOC01-appb-M000017
　上記の数式（１６）と数式（１７）とから、二乗平均平方根の二乗ｒｊ＋１ ２は、以下の数式（１８）のように表せる。
Figure JPOXMLDOC01-appb-M000018
　上記の数式（１８）より、第ｊ＋１時点の二乗平均平方根の二乗ｒｊ＋１ ２は、以下の数式（１９）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000019
　上記のように、第ｊ＋１時点の変数Ｌｊ＋１は、前述の数式（７）で逐次的に求められる。また、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第ｊ＋１時点の二乗平均平方根の二乗ｒｊ ２を求めるために第ｊ時点で保持すべき内部変数は、二乗平均平方根の二乗ｒｊ ２及び変数Ｌｊである。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ及び二乗平均平方根の二乗ｒｊ ２の初期値Ｌ１，ｒ１ ２は０とするか、初期化処理直後の変動を避けるため初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、特徴量の一例として、歪度ｗの逐次計算に関する計算手順を示す。センサデータｘを確率変数と見なし、確率変数ｘの平均値ｍ及び標準偏差ｓを用いると、歪度ｗは、以下の公知の数式（２０）で定義される。
Figure JPOXMLDOC01-appb-M000020
　前述したように、記号“Ｅ［］”は、角括弧内の確率変数の期待値を表している。上記の数式（２０）を展開し、整理すると、以下の数式（２１）のように表せる。
Figure JPOXMLDOC01-appb-M000021
　上記の数式（２１）の変形では、センサデータｘの期待値Ｅ［ｘ］が、平均値ｍと等しいという性質を利用している。
　次に、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する歪度ｗｊについて考える。上記数式（２１）において、右辺の分子第１項の期待値Ｅ［ｘ３］を変数Ａｊで表し、右辺の分子第２項の期待値を変数Ｂｊで表すと、歪度ｗｊは、以下の数式（２２）のように表せる。
Figure JPOXMLDOC01-appb-M000022
　上記の数式（２２）において、右辺の分母にある標準偏差ｓｊに関しては、前述の数式（１５）で逐次的に計算できる。
　次に、上記の数式（２２）の第ｊ時点の変数Ａｊ，Ｂｊに関する逐次計算を考える。まず、第ｊ時点の変数Ａｊは、以下の数式（２３）で表せる。
Figure JPOXMLDOC01-appb-M000023
　同様に、第ｊ＋１時点の変数Ａｊ＋１は、以下の数式（２４）で表せる。
Figure JPOXMLDOC01-appb-M000024
　上記の数式（２３）と数式（２４）とから、変数Ａｊ＋１は、以下の数式（２５）のように表せる。
Figure JPOXMLDOC01-appb-M000025
　従って、上記の数式（２５）より、第ｊ＋１時点の変数Ａｊ＋１は、以下の数式（２６）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000026
　上記のように、第ｊ＋１時点の変数Ｌｊ＋１は、前述の数式（７）で逐次的に求められる。また、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。ここで、第ｊ＋１時点の変数Ａｊ＋１を求めるために第ｊ時点で保持すべき内部変数は、変数Ａｊ及び変数Ｌｊである。
　また、第ｊ時点の変数Ｂｊは、以下の数式（２７）で表せる。
Figure JPOXMLDOC01-appb-M000027
　上記の数式（２７）のように、変数Ｂｊは前述の二乗平均平方根の二乗ｒｊ ２と一致することが分かり、これは前述の数式（１９）で逐次的に求められる。従って、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する歪度ｗｊは、以下の数式（２８）で求まる。
Figure JPOXMLDOC01-appb-M000028
　同様に、第１時点から第ｊ＋１時点までのセンサデータｘｊ＋１に対する歪度ｗｊ＋１は、以下の数式（２９）で求まる。
Figure JPOXMLDOC01-appb-M000029
　上記のように、第ｊ＋１時点の変数Ａｊ＋１は、前述の数式（２６）で逐次的に求められる。また、第ｊ＋１時点の平均値ｍｊ＋１は、前述の数式（６）で逐次的に求められる。そして、第ｊ＋１時点の二乗平均平方根の二乗ｒｊ ２は、前述の数式（１９）で逐次的に求められる。また、第ｊ＋１時点の標準偏差ｓｊ＋１は、前述の数式（１５）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第１時点から第ｊ＋１時点までのセンサデータｘｊ＋１に対する歪度ｗｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、変数Ａｊ，Ｌｊ、平均値ｍｊ、二乗平均平方根の二乗ｒｊ ２及び標準偏差ｓｊである。標準偏差ｓｊに代えて分散ｖｊを保持してもよい。また、二乗平均平方根の二乗ｒｊ ２に代えて二乗平均平方根ｒｊを保持してもよい。
　なお、数式が煩雑になるため割愛するが、逐次計算のための変数Ａｊに代えて、歪度ｗｊを保持しても同様の計算結果が得られる数式が求まる。また、第２時点の内部変数を計算する際、簡単のため、変数Ａｊ，Ｌｊ、及び平均値ｍｊ、二乗平均平方根の二乗ｒｊ ２、標準偏差ｓｊの各初期値は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、特徴量の一例として、尖度ｋの逐次計算に関する計算手順を示す。センサデータｘを確率変数と見なし、確率変数ｘの標準偏差ｓｊを用いると、尖度ｋは、以下の公知の数式（３０）で定義される。
Figure JPOXMLDOC01-appb-M000030
　前述したように、記号“Ｅ［］”は、角括弧内の確率変数の期待値を表している。なお、本稿で説明する計算方法では、正規分布となる確率変数に対する尖度は３である。本稿と異なる定義として、正規分布の尖度が０となるよう、数式（３０）の右辺に“－３”のバイアスをかける定義も存在するが、同様の手順で以下の数式に対応する数式の導出が可能である。
　上記の数式（３０）を展開し、整理すると、以下の数式（３１）のように表せる。
Figure JPOXMLDOC01-appb-M000031
　上記の数式（３１）の変形では、センサデータｘの期待値Ｅ［ｘ］が、平均値ｍと等しいという性質を利用している。
　次に、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する尖度ｋｊについて考える。上記数式（３１）において、右辺の分子第１項の期待値Ｅ［ｘ４］を変数Ｃｊで表す。また、右辺の分子第２項の期待値及び第３項の期待値は、それぞれ変数Ａｊ，Ｂｊで表せる。従って、尖度ｋｊは、以下の数式（３２）で表せる。
Figure JPOXMLDOC01-appb-M000032
　上記の数式（３２）において、標準偏差ｓｊ、変数Ａｊ，Ｂｊに関しては、それぞれ前述の数式（１５）、数式（２６）及び数式（２７）で逐次的に計算できる。
　次に、上記の数式（３２）の第ｊ時点の変数Ｃｊに関する逐次計算を考える。まず、第ｊ時点の変数Ｃｊは、以下の数式（３３）で表せる。
Figure JPOXMLDOC01-appb-M000033
　同様に、第ｊ＋１時点の変数Ｃｊ＋１は、以下の数式（３４）で表せる。
Figure JPOXMLDOC01-appb-M000034
　上記の数式（３３）と数式（３４）とから、変数Ｃｊ＋１は、以下の数式（３５）のように表せる。
Figure JPOXMLDOC01-appb-M000035
　上記の数式（３５）より、第ｊ＋１時点の変数Ｃｊ＋１は、以下の数式（３６）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000036
　ここで、上記数式（２７）のように、変数Ｂｊは、前述の二乗平均平方根の二乗ｒｊ ２と一致する。従って、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する尖度ｋｊは、以下の数式（３７）で求まる。
Figure JPOXMLDOC01-appb-M000037
　同様に、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する尖度ｋｊ＋１は、以下の数式（３８）で求まる。
Figure JPOXMLDOC01-appb-M000038
　上記のように、第ｊ＋１時点の変数Ｃｊ＋１は、前述の数式（３６）で逐次的に求められる。また、第ｊ＋１時点の変数Ａｊ＋１は、前述の数式（２６）で逐次的に求められる。そして、第ｊ＋１時点の平均値ｍｊ＋１は、前述の数式（６）で逐次的に求められる。また、第ｊ＋１時点の二乗平均平方根の二乗ｒｊ ２は、前述の数式（１９）で逐次的に求められる。また、第ｊ＋１時点の標準偏差ｓｊ＋１は、前述の数式（１５）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第１時点から第ｊ＋１時点までのセンサデータｘｊ＋１に対する尖度ｋｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、変数Ａｊ，Ｌｊ，Ｃｊ、平均値ｍｊ、二乗平均平方根の二乗ｒｊ ２、標準偏差ｓｊである。標準偏差ｓｊに代えて分散ｖｊを保持してもよい。また、二乗平均平方根の二乗ｒｊ ２に代えて二乗平均平方根ｒｊを保持してもよい。
　なお、第２時点の内部変数を計算する際、簡単のため、変数Ａｊ，Ｌｊ，Ｃｊ、及び平均値ｍｊ、二乗平均平方根の二乗ｒｊ ２、標準偏差ｓｊの各初期値は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。また、センサデータが正規分布に近い性質を持つことが分かっている場合、尖度は３程度の値になることが知られている。そのため、尖度ｋの初期値ｋ１を３としてもよい。
　次に、特徴量の一例として、最大値の逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する最大値ａｊは、以下の公知の数式（３９）で定義される。
Figure JPOXMLDOC01-appb-M000039
　最大値ａの逐次的な計算方法は、以下の方法で容易に実現できることが知られている。第ｊ時点の最大値ａｊを保持しておき、第ｊ＋１時点のセンサデータｘｊと第ｊ時点の最大値ａｊとのうちで大きな方を、第ｊ＋１時点の最大値ａｊ＋１として保持する。即ち、第ｊ＋１時点の最大値ａｊ＋１は、以下の数式（４０）で逐次的に計算できる。
Figure JPOXMLDOC01-appb-M000040
　なお、最大値の初期値、即ち最大値ａ１には、同時刻のセンサデータｘ１を代入するのがよい。
　次に、特徴量の一例として、最小値の逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する最小値ｎｊは、以下の公知の数式（４１）で定義される。
Figure JPOXMLDOC01-appb-M000041
　最小値ｎの逐次的な計算方法は、以下の方法で容易に実現できることが知られている。第ｊ時点の最小値ｎｊを保持しておき、第ｊ＋１時点のセンサデータｘｊと第ｊ時点の最小値ｎｊとのうちで小さな方を、第ｊ＋１時点の最小値ｎｊ＋１として保持する。即ち、第ｊ＋１時点の最小値ｎｊ＋１は、以下の数式（４２）で逐次的に計算できる。
Figure JPOXMLDOC01-appb-M000042
　なお、最小値の初期値、即ち最小値ｎ１には、同時刻のセンサデータｘ１を代入するのがよい。
　次に、特徴量の一例として、ピーク値の逐次計算に関する計算手順を示す。ピーク値という言葉には幾つかの定義が存在するが、本稿では、第１時点から第ｊ時点までのピーク値ｐｊを以下の数式（４３）で定義する。
Figure JPOXMLDOC01-appb-M000043
　上記数式（４３）の定義では、ピーク値はセンサデータｘ１～ｘｊのうちで絶対値の最も大きな値としている。ピーク値の逐次的な計算方法は、以下の方法で実現できる。第ｊ時点のピーク値ｐｊを保持しておき、第ｊ＋１時点のセンサデータｘｊと第ｊ時点のピーク値ｐｊのうちで大きな方を、第ｊ＋１時点のピーク値ｐｊ＋１として保持する。即ち、第ｊ＋１時点のピーク値ｐｊ＋１は、以下の数式（４４）で逐次的に計算できる。
Figure JPOXMLDOC01-appb-M000044
　ピーク値の初期値、即ちピーク値ｐ１には、同時刻のセンサデータの絶対値|ｘ１|を代入するのがよい。或いは、前述の数式（４０）及び数式（４２）を用いて、第ｊ時点の最大値ａｊと最小値ｎｊとが得られている場合、以下の数式（４５）でピーク値を求めることもできる。
Figure JPOXMLDOC01-appb-M000045
　同様に、第ｊ＋１時点の最大値ａｊ＋１と最小値ｎｊ＋１とが得られている場合、以下の数式（４６）でピーク値を求めることができる。
Figure JPOXMLDOC01-appb-M000046
　結果として、第ｊ＋１時点のピーク値ｐｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、ピーク値ｐｊである。なお、ピーク値ｐｊに代えて、最大値ａｊ及び最小値ｎｊを保持してもよい。
　次に、特徴量の一例として、ピークピーク値の逐次計算に関する計算手順を示す。ピークピーク値は、最大値と最小値との差である。ピークトゥピーク（ｐｅａｋ－ｔｏ－ｐｅａｋ）と呼ばれることもある。
　第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対するピークピーク値ｐｐｊは、以下の公知の数式（４７）で定義される。
Figure JPOXMLDOC01-appb-M000047
　即ち、前述の数式（４０）と数式（４２）とから、最大値ａｊ及び最小値ｎｊを逐次的に計算しておくことで、ピークピーク値ｐｐｊが直ちに求まる。
　同様に、第ｊ＋１時点の最大値ａｊ＋１と最小値ｎｊ＋１とが得られている場合、以下の公知の数式（４８）でピークピーク値ｐｐｊ＋１を求めることができる。
Figure JPOXMLDOC01-appb-M000048
　結果として、第ｊ＋１時点のピークピーク値ｐｐｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、最大値ａｊ及び最小値ｎｊである。
　次に、特徴量の一例として、波高率の逐次計算に関する計算手順を示す。第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する波高率ｐｒｊは、以下の数式（４９）で定義される。
Figure JPOXMLDOC01-appb-M000049
　上記数式（４９）に示されるように、波高率ｐｒｊは、前述のピーク値ｐｊを前述の二乗平均平方根ｒｊで割った値である。波高率ｐｒｊは、ピークトゥアールエムエス（ｐｅａｋ－ｔｏ－ＲＭＳ）、クレストファクタ（ｃｒｅｓｔ　ｆａｃｔｏｒ）と呼ばれることもある。なお、センサデータが正の方向に偏った値を取ることがわかっている場合などには、ピーク値ｐｊに代えて最大値ａｊを用いることもある。
　同様に、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊに対する波高率ｐｒｊ＋１は、以下の数式（５０）で求まる。
Figure JPOXMLDOC01-appb-M000050
　即ち、前述の数式（４４）又は数式（４５）を用いてピーク値ｐｊ＋１を逐次的に計算すると共に、数式（１９）を用いて二乗平均平方根ｒｊ＋１を逐次的に計算し逐次的に計算しておくことで、波高率ｐｒｊ＋１が直ちに求まる。
　結果として、第ｊ＋１時点の波高率ｐｒｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、ピーク値ｐｊ及び二乗平均平方根の二乗ｒｊ ２である。なお、ピーク値ｐｊに代えて、最大値ａｊ及び最小値ｎｊを保持してもよい。また、二乗平均平方根の二乗ｒｊ ２に代えて、二乗平均平方根ｒｊを保持してもよい。
　上述のように、内部変数計算部１００３は、各特徴量を計算するための内部変数を逐次計算する。内部変数計算部１００３は、この逐次計算をｊ＝１からｊ＝Ｎ―１まで実施する。これにより、特徴量計算部１００４が第１時点から第Ｎ時点でのセンサデータｘ１～ｘｊに対応する特徴量を決定することができる。なお、特徴量の種類毎に更新の数式、及び内部変数保持部１００２が保持する内部変数は異なる。その一方で、内部変数保持部１００２が内部変数を保持し、内部変数計算部１００３が内部変数を逐次計算し、特徴量計算部１００４が内部変数に基づいて特徴量を計算する点は、特徴量の種類に依らず同様である。
　次に、前述の方法を用いてセンサデータの一つであるモータトルクの特徴量を逐次的に計算した結果について、図６を参照して説明する。図６は、実施の形態１におけるモータトルク及び各種の特徴量の時系列の波形を示す図である。なお、以降では、記載の煩雑さを避けるため、特徴量等に付していた記号の表記を適宜省略する。
　図６（ａ）には、図５（ｂ）に示したものと同様のモータトルクの波形が示されている。図６（ａ）、（ｂ）の横軸は、時間を表している。また、図６（ａ）には、有次元の特徴量のうちモータトルクと同じ単位［Ｎｍ］を持つ、平均値、標準偏差、二乗平均平方根（ＲＭＳ）、最大値及び最小値の時系列の波形が示されている。またピークピーク値は、最大値と最小値との差であり、簡単のため図示しない。また、ピーク値は、最大値と最小値のうち、絶対値が大きいものを示す値であり、簡単のため図示しない。また、分散は、標準偏差を二乗したものであり、単位が異なるので図示しない。
　図６（ｂ）には、無次元の特徴量、即ち単位を持たない特徴量である、歪度、尖度及び波高率の時系列の波形が示されている。次に、図６（ａ）、（ｂ）に示される各特徴量の挙動に関して説明する。
　平均値は、期間Ｔｓ１ではほぼ０であり、期間Ｔａでは徐々に増加し、期間Ｔｅでは僅かに減少し、期間Ｔｄでは徐々に減少し、期間Ｔｓ２では僅かに減少している。
　標準偏差は、期間Ｔｓ１ではほぼ一定の小さい値であり、期間Ｔａでは徐々に増加し、期間Ｔｅでは僅かに減少し、期間Ｔｄでは徐々に増加し、期間Ｔｓ２では僅かに減少している。
　二乗平均平方根は、期間Ｔｓ１ではほぼ一定の小さい値であり、期間Ｔａでは徐々に増加し、期間Ｔｅでは僅かに減少し、期間Ｔｄではほぼ一定であり、期間Ｔｓ２では僅かに減少している。
　最大値は、期間Ｔｓ１ではほぼ一定の正の小さい値であり、期間Ｔａではモータトルクが過去の最大値を超える毎に増加しており、期間Ｔｅ、期間Ｔｄ及び期間Ｔｓ２では一定である。
　最小値は、期間Ｔｓ１ではほぼ一定の負の小さい値であり、期間Ｔａ及び期間Ｔｅでは一定、期間Ｔｄではモータトルクが過去の最小値を超える毎に減少しており、期間Ｔｓ２では一定である。
　歪度は、期間Ｔｓ１ではほぼ０であり、期間Ｔａが始まる時刻Ｔｒ１では急峻に大きい正の値に増加し、その後の期間Ｔａでは徐々に小さい値に減少し、期間Ｔｅでは僅かに増加し、期間Ｔｄでは僅かに減少し、期間Ｔｓ２では僅かに増加している。
　尖度は、期間Ｔｓ１では２程度から３程度まで増加し、期間Ｔａが始まる時刻Ｔｒ１では急峻に大きい正の値に増加し、その後の期間Ｔａでは徐々に小さい値に減少し、期間Ｔｅでは僅かに増加し、期間Ｔｄでは僅かに増加した後に減少し、期間Ｔｓ２では僅かに増加している。
　波高率は、期間Ｔｓ１では２程度から３程度まで増加し、期間Ｔａが始まる時刻Ｔｒ１では急峻に大きい正の値に増加し、その後の期間Ｔａでは徐々に小さい値に減少し、期間Ｔｅでは僅かに増加し、期間Ｔｄではほぼ一定であり、期間Ｔｓ２では僅かに増加している。
　図６（ａ）、（ｂ）に示した各時点の特徴量は、各時点より過去のモータトルクに対する特徴量が逐次的に計算されたものである。そのため、例えば平均値は、図６（ａ）の全期間（Ｔｓ１＋Ｔａ＋Ｔｅ＋Ｔｄ＋Ｔｓ２）のモータトルクを保持して一括計算で求めた値と、逐次的に求めた図６（ａ）の最終時刻である時刻Ｔｒ５の値とは一致する。
　なお、図６では、特徴量の推移を変わりやすくするために、特徴量計算部１００４もまた、センサデータの取得周期毎に各特徴量を計算している。ここで、ある期間に渡るセンサデータに対して特徴量を計算する場合、内部変数計算部１００３は、センサデータの取得周期毎に内部変数を逐次計算する必要がある。一方、例えば数式（２９）で表される歪度、数式（３８）で表される尖度のように、自身を内部変数として含めない特徴量を用いる場合、特徴量計算部１００４は、必要に応じて特徴量を計算することとしてもよい。例えば、時刻Ｔｒ０時点から時刻Ｔｒ５までの期間に渡るセンサデータに対する特徴量を１点だけ求めたく、且つ途中の時間での特徴量の推移が不要な場合である。このような場合、特徴量計算部１００４は、時刻Ｔｒ５の時点、即ちセンサデータのデータ点数Ｎが１２０００に達した時点の１回だけ特徴量を計算するように構成されていてもよい。このように構成すれば、情報処理装置１０００の計算負荷の低減が可能である。
　次に、状態診断部１００６が、特徴量を用いて機械装置１００８の状態を診断する手順について、図７を参照して説明する。図７は、実施の形態１における特徴量の経年変化の様子の一例を示す図である。図７には、前述の方法で計算した特徴量を定期的にプロットした結果が示されている。定期的とは、例えば一月毎である。図７の横軸は時間であり、機械装置１００８の運用時間を表している。
　状態診断部１００６は、特徴量に対して、しきい値Ｆｔｈ１を設け、特徴量がしきい値Ｆｔｈ１を超えた場合に、機械装置１００８に異常が生じていると診断する。特徴量には、前述したどのような特徴量を用いてもよい。また、複数の特徴量をそれぞれ計算し、各特徴量に個別のしきい値を設定してもよい。
　図７において、時刻Ｔｔａ０は運用を開始した時刻である。時刻Ｔｔａ０から時刻Ｔｔａ１より前までの特徴量は、多少のばらつきを生じながらも、しきい値Ｆｔｈ１よりも小さい値を保っている。状態診断部１００６は、特徴量がしきい値Ｆｔｈ１を下回っている場合には、機械装置１００８が正常である旨の診断を下し、その診断結果を出力する。
　また、図７の例では、時間の経過と共に特徴量の平均的な値が徐々に増加していき、時刻Ｔｔａ１で特徴量がしきい値Ｆｔｈ１を上回っている。この場合、状態診断部１００６は、時刻Ｔｔａ１で機械装置１００８に異常が生じたという診断を下し、その診断結果を出力する。
　しきい値Ｆｔｈ１の決め方は様々である。例えば過去の別の機械装置に異常が生じた際の特徴量に基づいて、しきい値Ｆｔｈ１を決定してもよい。また、機械装置１００８の運用を開始した直後の特徴量を参考にして、しきい値Ｆｔｈ１を決定してもよい。また、しきい値Ｆｔｈ１は、動的に設定してもよい。例えば複数の機械装置１００８を運用しながら、それらの特徴量を個別に計算し、装置ごとの特徴量のばらつきを考慮して定期的にしきい値Ｆｔｈ１を変更してもよい。また、機械装置１００８の特性を考慮したシミュレーション等で求めた特徴量に基づいて、しきい値Ｆｔｈ１を決定してもよい。また、本実施の形態において、しきい値Ｆｔｈ１は、正常な状態の上限値を表しているが、より小さな値が異常を示す特徴量を用いる場合は、正常な状態の下限値をしきい値Ｆｔｈ１として設定してもよい。更に、上限値及び下限値からなる２つのしきい値を併用してもよい。また、診断結果は一つである必要はなく、複数の特徴量に対して診断を実施し、特徴量毎に診断結果を出力してもよい。
　また、本実施の形態では、特徴量が一度でもしきい値Ｆｔｈ１を超えた場合に、機械装置１００８が異常であるという診断を下しているが、この手法に限定されない。ノイズ等の影響による誤った診断、即ち誤判定を減らすため、一定の期間の異常度の分布に基づいた検定（Ｓｔａｔｉｓｔｉｃａｌ　ｔｅｓｔ）を行うなどの統計的な手法を採用してもよい。
　以上に説明した処理により、図８に示すフローチャートが導かれる。図８は、実施の形態１における情報処理方法の説明に供するフローチャートである。
　まず、センサデータ取得部１００１は、機械装置１００８の物理量を計測するセンサ１０１０から、第１時点から第Ｎ時点（Ｎは２以上の整数）までの各時点における物理量の計測値をセンサデータとして取得する（ステップＳ１０１）。内部変数保持部１００２は、センサデータに基づいて時系列に順次計算される、Ｎよりも少ない個数の内部変数を保持する（ステップＳ１０２）。内部変数計算部１００３は、第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する内部変数を、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算する（ステップＳ１０３）。特徴量計算部１００４は、第Ｎ時点の内部変数に基づいて第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量を計算する（ステップＳ１０４）。状態診断部１００６は、特徴量に基づいて機械装置の状態を診断する（ステップＳ１０５）。
　なお、図８のフローチャートでは省略したが、第１時点の内部変数を予め設定した最大値と最小値との間の値に決定する初期化処理を実行するステップが含まれていてもよい。
　以上説明したように、実施の形態１における情報処理装置によれば、センサデータ取得部は、センサによって計測された機械装置の物理量の計測値を取得し、当該計測値のうちの第１時点から第Ｎ時点（Ｎは２以上の整数）までの計測値をセンサデータとして保持する。内部変数保持部は、センサデータに基づいて時系列に順次計算される、Ｎよりも少ない個数の内部変数を保持する。内部変数計算部は、第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する内部変数を、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算する。特徴量計算部は、第Ｎ時点の内部変数に基づいて第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量を計算する。このように構成された情報処理装置によれば、Ｎよりも少ない個数の内部変数を保持することにより、第１時点から第Ｎ時点までのセンサデータを保持する必要がない。これにより、計算能力の高い計算機及び記憶容量の大きな計算機を用いることなく、特徴量の計算処理が可能になるという従来にない顕著な効果が得られる。
　また、実施の形態１における情報処理装置は、特徴量に基づいて機械装置の状態を診断する状態診断部を備える。これにより、計算能力の高い計算機及び記憶容量の大きな計算機を用いることなく、機械装置の状態を診断するための情報処理を実施することができる。
　なお、実施の形態１における情報処理装置は、第１時点の内部変数を予め設定した最大値と予め設定した最小値との間の値に決定する初期化処理を実行する初期化処理部を備えていてもよい。このような初期化処理部を備えることにより、誤判定の確率を低減させ、診断処理の精度を高めることができる。
　また、実施の形態１における情報処理方法によれば、第１ステップでは、機械装置の物理量を計測するセンサから、第１時点から第Ｎ時点（Ｎは２以上の整数）までの各時点における物理量の計測値がセンサデータとして取得される。第２ステップでは、センサデータに基づいて時系列に順次計算される、Ｎよりも少ない個数の内部変数が保持される。第３ステップでは、第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する内部変数が、第ｊ＋１時点のセンサデータと第ｊ時点に対応する内部変数とに基づいて計算される。第４ステップでは、第Ｎ時点の内部変数に基づいて第１時点から第Ｎ時点までのセンサデータに含まれる統計的な特徴を抽出した特徴量が計算される。第５ステップでは、特徴量に基づいて機械装置の状態が診断される。このような、第１から第５ステップの処理を含む情報処理方法によれば、Ｎよりも少ない個数の内部変数が保持されるので、第１時点から第Ｎ時点までのセンサデータを保持する必要がない。これにより、計算能力の高い計算機及び記憶容量の大きな計算機を用いることなく、機械装置の状態を診断するための情報処理を実施できるという従来にない顕著な効果が得られる。
　なお、実施の形態１における情報処理方法には、上記第１から第５ステップの他、第１時点の内部変数を予め設定した最大値と予め設定した最小値との間の値に決定する初期化処理を実行するステップが含まれていてもよい。このような初期化処理のステップを含むことにより、誤判定の確率を低減させ、診断処理の精度を高めることができる。
実施の形態２．
　図９は、実施の形態２における情報処理装置２０００を含む情報処理システム１００Ａの構成例を示す図である。図９に示す情報処理システム１００Ａでは、図１に示す情報処理システム１００の構成において、情報処理装置１０００が情報処理装置２０００に置き替えられ、装置制御部１０９９が装置制御部２０９９に置き替えられている。情報処理装置２０００では、センサデータ取得部１００１がセンサデータ取得部２００１に置き替えられ、内部変数計算部１００３が内部変数計算部２００３に置き替えられ、特徴量計算部１００４が特徴量計算部２００４に置き替えられている。その他の構成は、図１に示す情報処理システム１００と同一又は同等である。なお、同一又は同等の構成部には同一の符号を付すと共に、重複する説明は省略する。
　センサデータ取得部２００１は、センサデータ取得部１００１の処理に加え、動作信号に基づいて内部変数を更新するか否かを決定するための計算許可フラグを生成する。より一般化して説明すると、センサデータ取得部２００１は、センサデータを取得する時間の間隔である取得周期、機械装置１００８の動作の状況を表す動作信号、センサデータにおけるデータ値のうちの少なくとも一つに基づいて第１時点から第ｊ時点（ｊは１からＮ－１までの整数（但し、Ｎは２以上の整数））までの間の２点以上の時点を決定し、決定した２点以上の時点に基づいて内部変数を更新するための計算許可フラグを生成する。
　内部変数計算部２００３は、内部変数計算部１００３の処理に加え、計算許可フラグに内部変数を更新する旨の内容が含まれている場合に内部変数を更新する。一方、計算許可フラグに内部変数を更新する旨の内容が含まれていない場合、内部変数計算部２００３は、内部変数の値に前回時刻のものを引き継がせる処理を行う。より一般化して説明すると、内部変数計算部２００３は、計算許可フラグの示す各第ｕ時点（ｕは１以上、且つｊ以下の整数（但し、ｊは１からＮ－１までの整数、且つＮは２以上の整数））について、第ｕ＋１時点の内部変数を、第ｕ時点の内部変数の値に決定する。
　特徴量計算部２００４は、特徴量計算部１００４の処理に加え、計算許可フラグに特徴量を更新する旨の内容が含まれている場合に特徴量を更新する。計算許可フラグに特徴量を更新する旨の内容が含まれていない場合、特徴量計算部２００４は、特徴量に前回時刻のものを引き継がせる処理を行う。
　情報処理装置２０００の外部にある装置制御部２０９９は、装置制御部１０９９の処理に加え、機械装置１００８の動作の状況を表す動作信号を出力する。機械装置１００８が、例えばモータにより駆動される機械装置の場合、速度が一定の場合の特徴量を用いたときにノイズが少ないセンサデータが得られ、機械装置１００８の状態を高精度に推定できることがある。そのため、モータ１００９の速度が一定となっている場合にのみ内部変数や特徴量の値を更新するように構成することで、機械装置１００８の状態診断に更に有用な特徴量を計算できるようになる。
　次に、前述の方法を用いてセンサデータの一つであるモータトルクの特徴量を逐次的に計算した結果について、図１０を参照して説明する。図１０は、実施の形態２におけるモータトルク及び各種の特徴量の時系列の波形を示す図である。
　図１０（ａ）には、図５（ａ）に示したモータ速度の時系列の波形と同じものが示されている。また、図１０（ｂ）には、センサデータ取得部２００１が決定した計算許可フラグの時系列の波形が示されている。
　計算許可フラグの結果の示し方の例として、内部変数を更新する場合にフラグ値を“１”とし、更新しない場合にフラグ値を“０”とすることが挙げられる。センサデータ取得部２００１は、速度が“０“ではなく、且つ一定である場合に内部変数が更新されるように、時刻Ｔｒ２から時刻Ｔｒ３までの期間Ｔｅで計算許可フラグのフラグ値を“１”としている。また、期間Ｔｅ以外、即ち期間Ｔｓ１、期間Ｔａ、期間Ｔｄ及び期間Ｔｓ２において、センサデータ取得部２００１は、計算許可フラグのフラグ値を“０”としている。なお、図９では、速度一定の例を５００［ｒ／ｍｉｎ］としている。
　速度が０の場合に内部変数を更新しない理由の一つに、モータ１００９が停止している際は機械装置１００８が動作せず、機械装置１００８の状態に応じた情報がセンサデータに現れないことが挙げられる。但し、機械装置１００８の構成によっては、モータ１００９が停止している際にも機械装置１００８の状態に応じた情報がセンサデータに現れる場合がある。このような場合には、停止時にも内部変数が更新されるよう構成されていてもよい。なお、計算許可フラグのフラグ値が“０”の際の動作は、内部変数及び特徴量が変化しなければよく、例えば今回の値に前回の値を上書きするように計算を実施しても、計算自体をスキップしてもよい。
　次に、計算許可フラグの決定方法として、機械装置１００８の動作の状況を表す動作信号から決定した例を示す。本例では、モータの速度及び加速度の２つを動作信号として用いる。モータの加速度の絶対値が小さく、且つ速度の絶対値が大きい場合、計算許可フラグのフラグ値を“１”にする。基準となる加速度及び速度は、機械装置１００８又はモータ１００９の諸元等を参考に、予め定めておくのがよい。動作信号は、装置制御部２０９９から得られる。
　また、モータ１００９の動作が予め決められている場合には、モータ速度を参照しなくてもよい。例えば、モータ１００９が動き始めてから予め定めた時刻になった際に計算許可フラグのフラグ値を“１”にし、更に予め定めた別の時刻の経過後に計算許可フラグのフラグ値を“０”にするようにしてもよい。これらの処理を行う場合、センサデータを取得する時間間隔であるサンプリング周期を参照することができる。
　図１０（ｃ）には、図５（ｂ）に示したものと同様のモータトルクの時系列波形が示されている。また、図１０（ｃ）には、図６（ａ）に示したものと同じ特徴量の時系列波形が示されている。また、図１０（ｄ）には、図６（ｂ）に示したものと同じ特徴量の時系列波形が示されている。
　図１０（ｃ）、（ｄ）の何れの特徴量も、期間Ｔｅを除く期間Ｔｓ１、期間Ｔａ、期間Ｔｄ及び期間Ｔｓ２では値が変化していない。これは、図１０（ｃ）の計算許可フラグが示すように、期間Ｔｅを除く期間では内部変数及び特徴量の更新が実施されないよう、内部変数計算部２００３及び特徴量計算部２００４が構成されているためである。一方、期間Ｔｅでは、何れの特徴量も随時更新されており、期間Ｔｅの終了時点である時刻Ｔｒ３までに、ほぼ一定値に収束している。このような構成によれば、実施の形態１の効果に加え、機械装置１００８の状態を更に高精度に推定することができる。
　以上説明したように、実施の形態２における情報処理装置によれば、センサデータ取得部は、センサデータを取得する時間の間隔である取得周期、機械装置の動作の状況を表す動作信号、及びセンサデータにおけるデータ値のうちの少なくとも一つに基づいて第１時点から第ｊ時点までの間の２点以上の時点を決定し、決定した２点以上の時点に基づいて計算許可フラグを生成する。内部変数計算部は、計算許可フラグの示す各第ｕ時点（ｕは１以上、且つ前記ｊ以下の整数）について、第ｕ＋１時点の内部変数を、第ｕ時点の内部変数の値に決定する。このように構成された情報処理装置によれば、機械装置の動作状況に応じてノイズが少ないセンサデータを得ることができる。これにより、実施の形態１の効果を得た上で、更に機械装置の状態を高精度に推定できるという効果が得られる。
実施の形態３．
　図１１は、実施の形態３における情報処理装置３０００を含む情報処理システム１００Ｂの構成例を示す図である。図１１に示す情報処理システム１００Ｂでは、図９に示す情報処理システム１００Ａの構成において、情報処理装置２０００が情報処理装置３０００に置き替えられている。情報処理装置３０００では、センサデータ取得部２００１がセンサデータ取得部３００１に置き替えられ、初期化処理部１００５が初期化処理部３００５に置き替えられている。その他の構成は、図９に示す情報処理システム１００Ａと同一又は同等である。なお、同一又は同等の構成部には同一の符号を付すと共に、重複する説明は省略する。
　センサデータ取得部３００１は、図１に示すセンサデータ取得部１００１の処理に加え、内部変数を初期化する時点を決めるための初期化トリガを生成する。より一般化し、且つより具体的に説明すると、センサデータ取得部３００１は、センサデータを取得する時間の間隔である取得周期、機械装置１００８の動作の状況を表す動作信号、及びセンサデータのうちの少なくとも一つに基づいて第１時点から第Ｎ時点までの間の１点以上の時点を決定し、決定した１点以上の時点に基づいて内部変数を初期化するための初期化トリガを生成する。
　初期化処理部３００５は、初期化処理部１００５の処理に加え、初期化トリガの示す時点に基づいて内部変数の初期化処理を実施する。
　機械装置１００８が、例えばモータ１００９により駆動される機械装置の場合、加速中、一定速度で回転中、減速中、停止中などといった運転状況毎に特徴量の傾向が異なる場合がある。このような場合、運転状況毎に特徴量の計算を切り分けた方が、機械装置１００８の状態を高精度に推定できることがある。そのため、モータ１００９の速度変化が変更される時点など、運転状況が変わる時点で内部変数を初期化することで、機械装置１００８の状態推定に更に有用な特徴量を計算できるようになる。
　次に、前述の方法を用いてセンサデータの一つであるモータトルクの特徴量を逐次的に計算した結果について、図１２を参照して説明する。図１２は、実施の形態３におけるモータトルク及び各種の特徴量の時系列の波形を示す図である。
　図１２（ａ）には、図５の（ａ）に示したモータ速度の時系列の波形と同じものが示されている。また、図１２（ｂ）には、センサデータ取得部３００１が決定した、初期化トリガの時系列の波形が示されている。
　初期化トリガの結果の示し方の例として、内部変数を初期化する場合に信号レベルを“１”とし、初期化しない場合に信号レベルを“０”とすることが挙げられる。センサデータ取得部３００１は、モータ速度の変化、即ち加速度が生じた時点で内部変数を初期化するように、時刻Ｔｒ１、時刻Ｔｒ２、時刻Ｔｒ３及び時刻Ｔｒ４で初期化トリガの信号レベルを“１”とし、それ以外の時刻で信号レベルを“０”としている。
　次に、初期化トリガの決定方法として、機械装置１００８の動作の状況を表す動作信号から決定した例を示す。本例では、モータのジャーク、即ち加速度の時間微分値を動作信号として用いる。モータのジャークの絶対値が大きい場合、初期化トリガの信号レベルを“１”にしている。基準となるジャークの大きさは、機械装置１００８又はモータ１００９の諸元等を参考に、予め定めておくのがよい。動作信号は、装置制御部２０９９から得られる。
　また、モータ１００９の動作が予め決められている場合には、モータ速度を参照せずともよい。例えば、モータ１００９が動き始めてから予め定めた時刻になった際に、瞬間的に初期化トリガの信号レベルを“１”とするように構成してもよい。或いは、特定の時間が経過する毎に初期化トリガの信号レベルを瞬間的に“１”にするよう構成してもよい。これらの処理を行う場合、センサデータを取得する時間間隔であるサンプリング周期を参照することができる。
　また、モータ１００９の動作開始から停止までの一連のセンサデータ毎に特徴量を計算したい場合、モータ１００９が動き出した時点で初期化処理を実施するように、初期化トリガを決定してもよい。このように構成することで、モータの移動距離が長い動作と、短い動作とを一律に１回の位置決め動作として扱うことができる。
　図１２（ｃ）には、図５（ｂ）に示したものと同様のモータトルクの時系列波形が示されている。また、図１２（ｃ）には、図６（ａ）に示したものと同じ特徴量の時系列波形が示されている。また、図１２（ｄ）には、図６（ｂ）に示したものと同じ特徴量の時系列波形が示されている。
　図１０（ｃ）、（ｄ）の何れの特徴量も、時刻Ｔｒ０、時刻Ｔｒ１、時刻Ｔｒ２、時刻Ｔｒ３及び時刻Ｔｒ４の直後に値が変化し始め、各期間Ｔｓ１、期間Ｔａ、期間Ｔｅ、期間Ｔｄ及び期間Ｔｓ２が終わるまでに、各特徴量はほぼ一定値に収束している。このような構成によれば、実施の形態１の効果に加え、機械装置１００８の状態を更に高精度に推定することができる。また、実施の形態２の効果に加え、モータ１００９が一定速度を持たない条件で運転する場合にも適用できるので、適用範囲を拡大することができる。
　以上説明したように、実施の形態３における情報処理装置によれば、センサデータ取得部は、センサデータを取得する時間の間隔である取得周期、機械装置の動作の状況を表す動作信号、及びセンサデータにおけるデータ値のうちの少なくとも一つに基づいて第１時点から第ｊ時点までの間の１点以上の時点を決定し、決定した１点以上の時点に基づいて内部変数を初期化するための初期化トリガを生成する。初期化処理部は、初期化トリガに基づいて内部変数を初期化する処理を実行する。このように構成された情報処理装置によれば、運転状況毎に特徴量の計算を切り分けることができる。これにより、実施の形態１の効果を得た上で、更に、機械装置の状態を高精度に推定できるという効果が得られる。また、モータが一定速度を持たない条件で運転する場合にも適用できるので、実施の形態２の効果を得た上で、更に運転条件に関する適用範囲を拡大できるという効果が得られる。
実施の形態４．
　図１３は、実施の形態４における情報処理装置４０００を含む情報処理システム１００Ｃの構成例を示す図である。図１３に示す情報処理システム１００Ｃでは、図１１に示す情報処理システム１００Ｂの構成において、情報処理装置３０００が情報処理装置４０００に置き替えられている。情報処理装置４０００では、内部変数計算部１００３が内部変数計算部４００３に置き替えられ、特徴量計算部１００４が特徴量計算部４００４に置き替えられ、状態診断部１００６が状態診断部４００６に置き替えられている。また、新たに収束度計算部４００７が設けられている。その他の構成は、図１１に示す情報処理システム１００Ｂと同一又は同等である。なお、同一又は同等の構成部には同一の符号を付すと共に、重複する説明は省略する。
　内部変数計算部４００３は、内部変数計算部１００３の処理に加え、忘却係数に基づいて内部変数を計算する。特徴量計算部４００４も同様に、特徴量計算部１００４の処理に加え、当該忘却係数に基づいて特徴量を計算する。忘却係数は、時刻が古いセンサデータが特徴量に与える影響よりも、時刻が新しいセンサデータが特徴量に与える影響の方が大きくなるように重みをつけるための係数であり、０より大きく１より小さい値を有している。
　収束度計算部４００７は、内部変数に基づいて、特徴量の計算の収束の度合いを定量的に示す指標である収束度を計算する。状態診断部４００６は、収束度計算部４００７が計算した収束度が一定の条件を満たす時点の特徴量を基に、機械装置１００８の状態を診断する。
　次に、内部変数計算部４００３及び特徴量計算部４００４の処理について詳細に説明する。特徴量の具体的な例として、指数移動平均値、指数移動分散、指数移動標準偏差、指数移動二乗平均平方根、指数移動歪度及び指数移動尖度を用いたときの、各特徴量の逐次計算方法を、その導出手順と共に示す。
　まず、特徴量の一例として、指数移動平均の逐次計算に関する計算手順を示す。前述した平均値とは異なり、時刻毎に持たせた重みを考慮した重み付き平均値が知られている。重み付き平均値の例として、忘却係数λを乗じた重みが付与された重み付きのセンサデータを平均する方法がある。忘却係数λは、前述したように、時刻が古いセンサデータに与える影響よりも、時刻が新しいセンサデータに与える影響の方が大きくなるように重み付けするための係数である。以降、この処理を指数移動平均と呼び、その値を指数移動平均値と呼ぶ。指数移動平均値は、指数平滑移動平均値と呼ばれることもある。また、“１”から忘却係数λを減じた“１－λ”を平滑化係数と呼ぶ場合もある。
　まず、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動平均値ｍ'ｊは、以下の数式（５１）で定義される。
Figure JPOXMLDOC01-appb-M000051
　例えば、λ＝０．９と設定した場合、第３時点の指数移動平均値ｍ'３は、以下の数式（５２）のように表せる。
Figure JPOXMLDOC01-appb-M000052
　忘却係数λを用いた指数移動平均では、上記のように古いデータよりも最近のデータにおける係数が大きくなっている（０．８１＜０．９＜１）。これは平均を算出する際に最近のデータ、即ち時刻の新しいデータを重視しているといえる。忘却係数λが０に近いほど古いデータを忘れ易くなり、逆に１に近いほど古いデータを忘れにくくなる。
　次に、指数移動平均値の逐次計算方法を考える。第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する指数移動平均値ｍ'ｊ＋１は、以下の数式（５３）のように変形できる。
Figure JPOXMLDOC01-appb-M000053
　即ち、第ｊ＋１時点の指数移動平均値ｍ'ｊ＋１は、以下の数式（５４）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000054
　ここで、センサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。また、第ｊ＋１時点の変数Ｌ'ｊは、以下の数式（５５）で表せる。
Figure JPOXMLDOC01-appb-M000055
　更に、変数Ｌ'ｊの逐次計算の方法を考える。第ｊ＋１時点の変数Ｌ'ｊ＋１は、以下の数式（５６）のように表せる。
Figure JPOXMLDOC01-appb-M000056
　即ち、第ｊ＋１時点の変数Ｌ'ｊ＋１は、以下の数式（５７）で逐次的に計算できる。
Figure JPOXMLDOC01-appb-M000057
　結果として、第ｊ＋１時点の指数移動平均値ｍ'ｊ＋１を計算するのに保持しておくべき内部変数は、第ｊ時点の指数移動平均値ｍ'ｊ及び変数Ｌ'ｊである。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ'ｊ及び指数移動平均値ｍ'ｊの初期値Ｌ'１，ｍ'１は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、予め設定する忘却係数λの決定方法に関して説明する。忘却係数λは０より大きく、１より小さい値に設定するのが望ましい。忘却係数λは、１に近い場合に過去の情報を忘れにくく、０に近い場合に過去の情報を忘れ易くなる。忘却係数λを１にすると、センサデータの情報を忘却せず、逐次計算の数式は実施の形態１と一致する。また、忘却係数を０にすると、時刻が一つ経過した際に全てのセンサデータの情報を忘れるため、２つ以前過去のセンサデータに対する特徴量は求まらない。多くの場合、忘却係数は１に近く、且つ１未満の値、例えば０．９から０．９９９程度の値に設定するのがよい。忘却係数λを用いた指数移動平均等を計算する場合、センサデータを特徴量に反映する重みの多くは、時定数と呼ぶ直近の短い期間に集中する。ここで時定数τは、以下の数式（５８）で計算できる。
Figure JPOXMLDOC01-appb-M000058
　上記数式（５８）において、「ｄｔ」はサンプリング周期を表している。また、直近の特定の時間帯のセンサデータに対して多くの重みを付けて特徴量を計算したい場合には、「λ＝１－ｄｔ／τ」の数式に用いて忘却係数λを決定してもよい。この数式は、上記数式（５８）を変形したものである。
　例えば、直近の２０ミリ秒のセンサデータに対して多くの重みを付けて特徴量を計算したい場合には、時定数を２０ｍｓとする。このとき、サンプリング周期が０．５ｍｓであれば、忘却係数λ＝０．９７５となる。また、時定数τに代えて、時定数τの逆数である推定周波数（単位ｒａｄ／ｓ）を用いて忘却係数λを設定することとしてもよい。
　次に、特徴量の一例として、指数移動分散の逐次計算に関する計算手順を示す。前述のように、時刻毎に忘却係数λを乗じて重みを付加された重み付きセンサデータを用いた分散を「指数移動分散」と呼ぶ。
　まず、センサデータを確率変数と見なすと、その分散は前述の数式（９）で表せるということが知られている。そのため、第ｊ時点の指数移動分散ｖ'ｊは、以下の数式（５９）で表せる。
Figure JPOXMLDOC01-appb-M000059
　従って、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する指数移動分散ｖ'ｊ＋１は、以下の数式（６０）で表せる。
Figure JPOXMLDOC01-appb-M000060
　上記の数式（６０）を変形すると、指数移動分散ｖ'ｊ＋１は、以下の数式（６１）で表せる。
Figure JPOXMLDOC01-appb-M000061
　上記の数式（５９）と数式（６１）とから、以下の数式（６２）が導かれる。
Figure JPOXMLDOC01-appb-M000062
　従って、上記の数式（６２）より、第ｊ＋１時点の指数移動分散ｖ'ｊ＋１は、以下の数式（６３）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000063
　上記のように、第ｊ＋１時点の変数Ｌ'ｊ＋１は、前述の数式（５７）で逐次的に求められる。また、第ｊ＋１時点の指数移動平均値ｍ'ｊ＋１は、前述の数式（５４）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第ｊ＋１時点の指数移動分散ｖ'ｊ＋１を求めるために第ｊ時点で保持すべき内部変数は、変数Ｌ'ｊ、指数移動分散ｖ'ｊ及び指数移動平均値ｍ'ｊである。なお、指数移動分散ｖ'ｊに代えて、後述の指数移動標準偏差ｓ'ｊを保持してもよい。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ'、指数移動分散ｖ'及び指数移動平均ｍ'の初期値Ｌ'１，ｖ'１，ｍ'１は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、特徴量の一例として、指数移動標準偏差ｓ'ｊの逐次計算に関する計算手順を示す。前述のように、時刻毎に忘却係数λを乗じて重みを付加された重み付きセンサデータを用いた標準偏差を「指数移動標準偏差」と呼ぶ。
　まず、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動標準偏差ｓ'ｊは、指数移動分散ｖ'ｊから以下の数式（６４）で定義される。
Figure JPOXMLDOC01-appb-M000064
　即ち、指数移動標準偏差ｓ'ｊは、前述の手順で逐次的に求めた指数移動分散ｖ'ｊから直ちに求まる。そのため、指数移動標準偏差ｓ'ｊを逐次的に計算するために保持すべき内部変数は、指数移動分散ｖ'ｊと同じである。
　次に、特徴量の一例として、指数移動二乗平均平方根の逐次計算に関する計算手順を示す。前述のように、時刻毎に忘却係数λを乗じて重みを付加された重み付きセンサデータを用いた二乗平均平方根を「指数移動二乗平均平方根」と呼ぶ。
　まず、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動二乗平均平方根ｒ'ｊは、以下の数式（６５）で定義される。
Figure JPOXMLDOC01-appb-M000065
　従って、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する指数移動二乗平均平方根ｒ'ｊ＋１は、以下の数式（６６）のように表せる。
Figure JPOXMLDOC01-appb-M000066
　上記の数式（６５）と数式（６６）とから、指数移動二乗平均平方根の二乗ｒ'ｊ＋１ ２は、以下の数式（６７）のように表せる。
Figure JPOXMLDOC01-appb-M000067
　上記の数式（６７）より、第ｊ＋１時点の指数移動二乗平均平方根の二乗ｒ'ｊ＋１ ２は、以下の数式（６８）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000068
　上記のように、第ｊ＋１時点の変数Ｌ'ｊ＋１は、前述の数式（５７）で逐次的に求められる。また、センサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第ｊ＋１時点の指数移動二乗平均平方根の二乗ｒ'ｊ＋１ ２を求めるために第ｊ時点で保持すべき内部変数は、指数移動二乗平均平方根の二乗ｒ'ｊ ２及び変数Ｌ'ｊである。
　第２時点の内部変数を計算する際、簡単のため、変数Ｌ'ｊ及び指数移動二乗平均平方根の二乗ｒ'ｊ ２の初期値Ｌ'１，ｒ'１は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、特徴量の一例として、指数移動歪度の逐次計算に関する計算手順を示す。前述のように、時刻毎に忘却係数λを乗じて重みを付加された重み付きセンサデータを用いた歪度を「指数移動歪度」と呼ぶ。
　前述の数式（２０）より求められた数式（２２）を参照すると、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動歪度ｗ'ｊは、以下の数式（６９）で表せる。
Figure JPOXMLDOC01-appb-M000069
　また、上記の数式（６９）の指数移動平均値ｍ'ｊに関しては前述の数式（５４）で逐次的に計算でき、上記の数式（６９）の指数移動標準偏差ｓ'ｊに関しては前述の数式（６３）、（６４）で逐次的に計算できる。
　次に上記数式（６９）の第ｊ時点の変数Ａ'ｊ，Ｂ'ｊに関する逐次計算を考える。まず、第ｊ時点の変数Ａ'ｊ，Ｂ'ｊは、以下の数式（７０）、（７１）で表せる。
Figure JPOXMLDOC01-appb-M000070
Figure JPOXMLDOC01-appb-M000071
　また、第ｊ＋１時点の変数Ａ'ｊ＋１は、以下の数式（７２）で表せる。
Figure JPOXMLDOC01-appb-M000072
　上記の数式（７２）から、変数Ａ'ｊ＋１は、以下の数式（７３）ように表せる。
Figure JPOXMLDOC01-appb-M000073
　上記の数式（７３）より、第ｊ＋１時点の変数Ａ'ｊ＋１は、以下の数式（７４）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000074
　ここで、第ｊ＋１時点の変数を求めるために第ｊ時点で保持すべき内部変数は、変数Ｌ'ｊ，Ａ'ｊである。また、第ｊ＋１時点の変数Ｌ'ｊ＋１は前述の数式（５７）で逐次的に求められる。また、第ｊ時点の変数Ｂ'ｊは、以下の数式（７５）で表せる。
Figure JPOXMLDOC01-appb-M000075
　ここで、上記の数式（７５）のように、変数Ｂ'ｊは、前述の指数移動二乗平均平方根の二乗ｒ'ｊ ２と一致することが分かり、これは前述の数式（６８）で逐次的に求められる。従って、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動歪度ｗ'ｊは、以下の数式（７６）で求まる。
Figure JPOXMLDOC01-appb-M000076
　同様に、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する指数移動歪度ｗ'ｊ＋１は、以下の数式（７７）で求まる。
Figure JPOXMLDOC01-appb-M000077
　上記のように、第ｊ＋１時点の変数Ａ'ｊ＋１は、前述の数式（７４）で逐次的に求められる。また、第ｊ＋１時点の指数移動平均値ｍ'ｊ＋１は前述の数式（５４）で逐次的に求められる。そして、第ｊ＋１時点の指数移動二乗平均平方根の二乗ｒ'ｊ ２は、前述の数式（６８）で逐次的に求められる。また、第ｊ＋１時点の指数移動標準偏差ｓ'ｊ＋１は前述の数式（６４）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第１時点から第ｊ＋１時点までのセンサデータｘｊ＋１に対する指数移動歪度ｗ'ｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、変数Ａ'ｊ，Ｌ'ｊ、指数移動平均値ｍ'ｊ、指数移動二乗平均平方根の二乗ｒ'ｊ ２及び指数移動標準偏差ｓ'ｊである。指数移動標準偏差ｓ'ｊに代えて指数移動分散ｖ'ｊを保持してもよい。また、指数移動二乗平均平方根の二乗ｒ'ｊ ２に代えて指数移動二乗平均平方根ｒ'ｊを保持してもよい。
　なお、数式が煩雑になるため割愛するが、逐次計算のための変数Ａ'ｊに代えて、指数移動歪度ｗ'ｊを保持しても同様の計算結果が得られる数式が求まる。また、第２時点の内部変数を計算する際、簡単のため、変数Ａ'ｊ，Ｌ'ｊ、及び指数移動平均値ｍ'ｊ、指数移動二乗平均平方根の二乗ｒ'ｊ ２、指数移動標準偏差ｓ'ｊの各初期値は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　次に、特徴量の一例として、指数移動尖度の逐次計算に関する計算手順を示す。前述のように、時刻毎に忘却係数λを乗じて重みを付加された重み付きセンサデータを用いた尖度を「指数移動尖度」と呼ぶ。
　前述の数式（３０）より求められた数式（３２）を参照すると、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動尖度ｋ'ｊは、以下の数式（７８）で表せる。
Figure JPOXMLDOC01-appb-M000078
　ここで、第ｊ時点の変数Ｃ'ｊは、以下の数式（７９）で表せる。
Figure JPOXMLDOC01-appb-M000079
　また、上記の数式（７８）の変数Ａ'ｊ，Ｂ'ｊ、指数移動平均値ｍ'ｊ、指数移動標準偏差ｓ'ｊに関しては、前述の数式（５４）、（６３）、（６４）、（６８）、（７４）、（７５）で逐次的に計算できる。
　次に上記数式（７９）の第ｊ時点の変数Ｃ'ｊに関する逐次計算を考える。まず、第ｊ＋１時点の変数Ｃ'ｊ＋１は、以下の数式（８０）で表せる。
Figure JPOXMLDOC01-appb-M000080
　上記の数式（７９）と数式（８０）とから、変数Ｃ'ｊ＋１は、以下の数式（８１）のように表せる。
Figure JPOXMLDOC01-appb-M000081
　上記の数式（８１）より、第ｊ＋１時点の変数Ｃ'ｊ＋１は、以下の数式（８２）で逐次的に求められる。
Figure JPOXMLDOC01-appb-M000082
　ここで、上記の数式（７５）のように、変数Ｂ'ｊは、前述の指数移動二乗平均平方根の二乗ｒ'ｊ ２と一致し、指数移動二乗平均平方根の二乗ｒ'ｊ ２は、前述の数式（６８）で逐次的に求められる。従って、第１時点から第ｊ時点までのセンサデータｘ１～ｘｊに対する指数移動尖度ｋ'ｊは、以下の数式（８３）で求まる。
Figure JPOXMLDOC01-appb-M000083
　同様に、第１時点から第ｊ＋１時点までのセンサデータｘ１～ｘｊ＋１に対する指数移動尖度ｋ'ｊ＋１は、以下の数式（８４）で求まる。
Figure JPOXMLDOC01-appb-M000084
　上記のように、第ｊ＋１時点の変数Ｃ'ｊ＋１は、前述の数式（８２）で逐次的に求められる。また、第ｊ＋１時点の変数Ａ'ｊ＋１は前述の数式（７４）で逐次的に求められる。そして、第ｊ＋１時点の指数移動平均値ｍ'ｊ＋１は前述の数式（５４）で逐次的に求められる。また、第ｊ＋１時点の指数移動二乗平均平方根の二乗ｒ'ｊ＋１ ２は前述の数式（６８）で逐次的に求められる。また、第ｊ＋１時点の指数移動標準偏差ｓ'ｊ＋１は前述の数式（６３）、（６４）で逐次的に求められる。そして、第ｊ＋１時点のセンサデータｘｊ＋１は、第ｊ＋１時点で新しく取得される。結果として、第１時点か指数移動二乗平均平方根の二乗ｒ'ｊ ２は、前述の数式（６８）で逐次的に求められる。ら第ｊ＋１時点までのセンサデータｘｊ＋１に対する指数移動尖度ｋ'ｊ＋１を求めるために、第ｊ時点で保持すべき内部変数は、変数Ａ'ｊ，Ｌ'ｊ、指数移動平均値ｍ'ｊ、指数移動二乗平均平方根の二乗ｒ'ｊ ２及び指数移動標準偏差ｓ'ｊである。指数移動標準偏差ｓ'ｊに代えて指数移動分散ｖ'ｊを保持してもよい。また、指数移動二乗平均平方根の二乗ｒ'ｊ ２に代えて指数移動二乗平均平方根ｒ'ｊを保持してもよい。
　なお、数式が煩雑になるため割愛するが、逐次計算のための変数Ａ'ｊに代えて、指数移動歪度ｗ'ｊを保持しても同様の計算結果が得られる数式が求まる。また、第２時点の内部変数を計算する際、簡単のため、変数Ａ'ｊ，Ｌ'ｊ、及び指数移動平均値ｍ'ｊ、指数移動二乗平均平方根の二乗ｒ'ｊ ２、指数移動標準偏差ｓ'ｊの各初期値は０とするか、初期化処理直後の変動を避けるため、初期化処理時のセンサデータの値等に応じて決定するのがよい。
　また、センサデータが正規分布に近い性質を持つことが分かっている場合、指数移動尖度ｋ'ｊは３程度の値になることが知られている。そのため、指数移動尖度ｋ'ｊの初期値ｋ'１を３としてもよい。
　また、直近の時間に重みをつけた最大値、最小値、ピークピーク値、ピーク値及び波高率（以下、適宜「最大値等」と呼ぶ）の計算に関しては、実施の形態１で示したような、前回の最大値等から今回の最大値等を更新することはできない。しかしながら、例えば最大値等となる可能性が高い候補値を内部変数として幾つか保持しておき、規定の時刻が経過したら保持した候補値を削除するという手順を踏めばよい。このようにすれば、直近の時間のセンサデータに対応する最大値等を計算することができる。
　次に、収束度計算部４００７の動作に関する計算手順を示す。忘却係数λを用いた特徴量の逐次計算では、内部変数から逐次計算の収束度合いを定量的に表す収束度を計算することができる。例えば、内部変数の一つである変数Ｌ'ｊは、数式（５７）で逐次的に計算されることから、入力されるセンサデータに関係なく、１／（１－λ）の値へ漸近することが分かる。そのため、変数Ｌ'ｊの初期値Ｌ'１を０以上、且つ１／（１－λ）未満とした場合に、収束度が０から１までの値をとれるように、以下の数式（８５）で収束度ｄｉを定義することが考えられる。
Figure JPOXMLDOC01-appb-M000085
　本実施の形態では、上記の数式（８５）を使って収束度を計算するが、収束度は上記に限るものではなく、計算の収束に伴い増加する値であればどのような定義のものでもよい。但し、計算の収束に伴って値が単調増加し、特定の値に漸近するものが基準として利用しやすい。例えば、上記の数式（８５）の右辺全体をＭ乗（Ｍは正の実数）した、｛（１－λ）Ｌ'ｊ｝Ｍを収束度と定義する。この収束度は、０から１までの値をとり、且つ計算の収束に伴って単調増加するので好ましい定義式である。
　次に、前述の方法を用いてセンサデータの一つであるモータトルクの特徴量を逐次的に計算した結果について、図１４を参照して説明する。図１４は、実施の形態４におけるモータトルク及び各種の特徴量の時系列の波形を示す図である。
　図１４（ａ）には、図５の（ａ）に示したモータ速度の時系列の波形と同じものが示されている。また、図１４（ｂ）には、図１２（ｂ）に示した初期化トリガの時系列の波形と同じものが示されている。
　図１４（ｃ）には、収束度計算部４００７が数式（８５）を用いて計算した収束度の時系列の波形が示されている。逐次計算が始まった時刻Ｔｒ０、初期化トリガが“１”を示した、時刻Ｔｒ１、時刻Ｔｒ２、時刻Ｔｒ３及び時刻Ｔｒ４で、収束度は０になっている。これらは、初期化処理が実行され、内部変数の一つである第ｊ時点の変数Ｌ'ｊが０になったためである。ここで、第ｊ時点は、初期化トリガが“１”を示した時点を指している。収束度は、各初期化処理が実施された以降に増加し、徐々に１に漸近している。
　図１４（ｄ）には、図５の（ｂ）に示したものと同様のモータトルクの時系列波形が示されている。また、図１４（ｄ）には、モータトルクと同じ単位［Ｎｍ］を持つ、指数移動平均値、指数移動標準偏差及び指数移動二乗平均平方根（ＲＭＳ）の時系列の波形を示している。なお、図が煩雑になるのを避けるため「指数移動」の表記は省略している。
　図１４（ｅ）には、無次元の特徴量、即ち単位を持たない特徴量である、指数移動歪度及び指数移動尖度の時系列の波形を示している。ここで、忘却係数λは０．９７５としている。また、図１４（ｄ）と同様に、「指数移動」の表記は省略している。次に、図１４（ｄ）、（ｅ）の逐次計算された各特徴量の挙動に関して説明する。
　指数移動平均値は、時刻Ｔｒ０から時刻Ｔｒ５までの全期間に渡って、モータトルクに僅かに遅れつつも、モータトルクに追従している。これは、直近の短時間のモータトルクに対して大きな重みを付けて平均値を計算する本実施の形態の性質を示している。
　指数移動標準偏差は、期間Ｔｓ１では、ほぼ一定の小さい値になっている。期間Ｔａ、期間Ｔｅ、期間Ｔｄ及び期間Ｔｓ２では、初期化トリガが１になった直後に増加し、その後、次の初期化トリガが１になるまでの間減少し、ほぼ一定値に収束している。
　指数移動二乗平均平方根（ＲＭＳ）は、期間Ｔｓ１ではほぼ一定の小さい値になっている。期間Ｔａでは、モータトルクの増加に伴い徐々に増加している。期間Ｔｅでは、ほぼ一定値に収束している。期間Ｔｄでは、モータ１００９の絶対値の増大に伴い徐々に増加している。期間Ｔｓ２では、ほぼ一定の小さい値に収束している。
　指数移動歪度は、期間Ｔｓ１ではほぼ０になっている。これは、期間Ｔｓ１のモータトルクが左右対称の確率分布となっているためである。期間Ｔａが始まる時刻Ｔｒ１では、歪度は一時的に正の値になり、その後負の値になっている。これは、モータトルクが０から増加したため、０付近の古いデータと増大後の新しいデータとが一時的に混在し、モータトルクの確率分布が乱れているためである。
　また、指数移動歪度は、期間Ｔｅでは時刻Ｔｒ２の直後に一時的に正の値になり、その後ほぼ０に収束している。指数移動歪度は、期間Ｔｄでは時刻Ｔｒ３の直後に一時的に正の値になり、その後ほぼ一定の値に収束している。指数移動歪度は、期間Ｔｓ２では時刻Ｔｒ４の直後に一時的に負の値になり、その後ほぼ０に収束している。
　更に、指数移動尖度は、期間Ｔｓ１ではほぼ３程度の値になっている。値が３になるのは、正規分布に対する指数移動尖度の特性であり、期間Ｔｓ１でモータトルクが正規分布に近い性質をもっていたことを表している。指数移動尖度は、期間Ｔａが始まる時刻Ｔｒ１の直後に負の値となり、その後の正の値になっている。更にその後、指数移動尖度は、時刻Ｔｒ２までにほぼ一定の値に収束している。これは、モータトルクが０から増加したため、０付近の古いデータと増大後の新しいデータとが一時的に混在し、モータトルクの確率分布が乱れているためである。
　また、指数移動尖度は、期間Ｔｅが始まる時刻Ｔｒ２の直後に負の値となり、その後に正の値になっている。更にその後、指数移動尖度は、時刻Ｔｒ３までにほぼ一定の値に収束している。
　また、指数移動尖度は、期間Ｔｄが始まる時刻Ｔｒ３の直後に負の値となり、その後に正の値になっている。更にその後、指数移動尖度は、時刻Ｔｒ４までにほぼ一定の値に収束している。
　また、指数移動尖度は、期間Ｔｓ２が始まる時刻Ｔｒ４の直後に負の値となり、その後に正の値になっている。更にその後、指数移動尖度は、時刻Ｔｒ５までに徐々に一定の値に収束しつつある。
　次に、収束度計算部４００７によって計算された収束度の挙動、及び状態診断部４００６の動作に関して、図１５を用いて説明する。図１５は、図１４に示す初期化トリガ及び収束度の時系列の波形を時刻Ｔｒ２から時刻Ｔｒ３までの期間に渡って拡大した図である。
　図１５において、時刻Ｔｒ２では初期化トリガが“１”を示しており、初期化処理が実行される。初期化処理が実行されると、内部変数の一つである第ｊ時点の変数Ｌ'ｊが０になり、収束度が０になる。ここで、第ｊ時点は初期化トリガが“１”を示した時点を指している。収束度は、初期化後に増加し、徐々に１に漸近している。
　状態診断部４００６は、収束度計算部４００７が数式（８５）を用いて計算した収束度が一定の条件を満たす時点での特徴量を基に、機械装置１００８の状態を診断する。例えば、収束度に対してしきい値Ｃｔｈを設定する。状態診断部４００６は、収束度がしきい値Ｃｔｈを超えている時点の特徴量を用いて機械装置１００８の状態を診断する。ここで、しきい値Ｃｔｈは、収束度が漸近する１に近い値に設定するのがよい。しきい値Ｃｔｈの一例は０．９９である。
　ここで、しきい値Ｃｔｈを超えた時点の時刻を「Ｔｘ」と表記する。また、収束度がしきい値Ｃｔｈを超えない時刻Ｔｒ２から時刻Ｔｘまでの期間を「Ｔｎ」と表記する。また、収束度がしきい値Ｃｔｈを超えた時刻ＴｘからＴｒ３までの期間を「Ｔｃ」と表記する。
　状態診断部４００６は、期間Ｔｎにおいて計算された特徴量は機械装置１００８の診断には用いず、期間Ｔｃにおいて計算された特徴量を機械装置１００８の診断に用いるのがよい。このような構成によれば、内部変数が初期化された場合に計算の結果が十分収束したときの特徴量を用いることができる。これにより、実施の形態３の効果に加え、機械装置の状態を更に高精度に推定できる。
　以上説明したように、実施の形態４における情報処理装置は、実施の形態３の構成に対し、更に収束度計算部を備える。収束度計算部は、内部変数に基づいて、特徴量の計算の収束の度合いを定量的に示す指標である収束度を計算する。この構成により、内部変数が初期化された場合に計算の結果が十分収束したときの特徴量を用いることができる。これにより、実施の形態３の効果を得た上で、更に機械装置の状態を高精度に推定できるという効果が得られる。
実施の形態５．
　図１６は、実施の形態５における情報処理装置５０００を含む情報処理システム１００Ｄの構成例を示す図である。図１６に示す情報処理システム１００Ｄでは、図１に示す情報処理システム１００の構成において、状態診断部１００６が状態診断部５００６に置き替えられている。その他の構成は、図１に示す情報処理システム１００と同一又は同等である。同一又は同等の構成部には同一の符号を付すと共に、重複する説明は省略する。なお、本実施の形態では、状態診断部５００６を図１に示す情報処理装置１０００に適用する例で説明するが、図９、図１１又は図１３に示す、情報処理装置２０００，３０００，４０００のうちの何れに適用してもよい。
　図１７は、実施の形態５における状態診断部５００６の構成例を示す図である。状態診断部５００６は、状態量取得部５４０１と、学習部５４０２と、異常度計算部５４０３と、意思決定部５４０４とを備える。状態量取得部５４０１は、特徴量を含む状態量を取得する。学習部５４０２は、機械装置１００８が正常な状態の状態量に基づいて、機械装置１００８の状態と特徴量との関係を学習する。異常度計算部５４０３は、学習部５４０２が学習した学習結果に基づいて機械装置１００８の異常度合いを定量的に示す指標である異常度を計算する。意思決定部５４０４は、異常度計算部５４０３が計算した異常度に基づいて機械装置１００８の状態を診断した診断結果を決定する。
　次に、状態量取得部５４０１、学習部５４０２、異常度計算部５４０３及び意思決定部５４０４の動作について説明する。
　状態量取得部５４０１は、特徴量と共に、機械装置１００８に付与する設定値を含む情報を状態量として取得する。特徴量の一例として、前述の方法で逐次計算した歪度を用いることが考えられる。また、状態量の一つとして用いる機械装置１００８の情報として、モータ速度の設定値を利用することが考えられる。図５（ａ）に示した例では、モータ速度の最大が５００［ｒ／ｍｉｎ］となるよう装置制御部１０９９を設定している。本実施の形態では、モータ速度の設定値は５００［ｒ／ｍｉｎ］であるとする。なお、本実施の形態では、状態量取得部５４０１への入力情報が、歪度及びモータ速度の設定値という２つの場合を例示するが、この例に限定されない。状態量取得部５４０１への入力情報は、３以上であってもよい。
　図１８は、実施の形態５における特徴量とモータ速度の設定値との関係を示す図である。図１８には、特徴量の一例である歪度とモータ速度の設定値との組を１つのデータとして複数プロットしている。横軸はモータ速度の設定値であり、縦軸は歪度である。特徴量は複数でもよく、状態量は３次元の概念で表示されるものでもよい。また、モータ速度以外の状態量を用いてもよい。黒塗りの丸は、機械装置１００８が正常な状態である正常データである。白塗りの丸は、機械装置１００８に何らかの異常が発生した場合の異常データである。図１８に示すように正常データとして、機械装置１００８が正常な状態の状態量を数多く得ていることが望ましい。図１８をみると、正常データは、モータ速度の設定値が大きいほど特徴量が大きくなる傾向にあることがわかる。一方、異常データは、モータ速度の設定値が比較的小さいにも関わらず、特徴量が比較的大きくなっている。
　学習部５４０２及び異常度計算部５４０３の構成例として、主成分分析（ｐｒｉｎｃｉｐａｌ　ｃｏｍｐｏｎｅｎｔ　ａｎａｌｙｓｉｓ）を用いた異常度の計算方法が知られている。主成分分析は、多次元のデータである状態量を分散の大きい方向から順に軸を取り直す方法である。主成分分析を用いた学習部５４０２は、予め得られている正常データの分散共分散行列の固有値及び固有ベクトルを計算し、元の状態量を主成分の空間に線形写像する。元の状態量をベクトルｘ、写像後のベクトルをｙとすると、ベクトルｘの写像は、以下の数式（８６）で表せる。
Figure JPOXMLDOC01-appb-M000086
　学習部５４０２は、正常データを基に、この線形写像を表す行列である表現行列Ａを決定する。学習部５４０２は、複数用意した正常データから求めた分散共分散行列に関する固有値及び固有ベクトルが複数得られる。得られた複数の固有値に関し、固有値が大きい順に対応する複数の固有ベクトルを並べる。そして、並べて作成された複数の固有ベクトルによる行列が表現行列Ａとなる。本実施の形態の例で説明すると、２つの状態量を用いて２行１列の固有ベクトルが作成される。そして２行１列の固有ベクトルが２つ並べられて、２行２列の表現行列が作成される。なお、本例とは異なるが、得られた固有値の数より少ない固有ベクトルを並べて、表現行列を作成することで、写像後のベクトルの次元を、元の状態量のベクトルの次元よりも小さくしてもよい。
　図１９は、実施の形態５における主成分分析の結果を示す図である。図１９には、図１８のデータが再掲され、且つ固有値が最も大きい第１主成分の軸と、固有値が２番目に大きい第２主成分の軸とが示されている。固有値は、主成分分析の結果として得られる。また、図１９には、主成分分析の結果を分かりやすくするために、第１主成分の軸と第２主成分の軸が交差する位置を中心に、各主成分の分散に対応する固有値に従って、楕円状の信頼区間が示されている。信頼区間は、データの分布がある確率で収まると期待される区間である。例えば、正常データの９０％が収まると期待される９０％信頼区間、正常データの９９％が収まると期待される９９％信頼区間などがある。信頼区間は、一般的な統計学の手法で計算できる。
　図１９において、正常データは９９％信頼区間の内部にあるのに対し、異常データは９９％信頼区間の外部に位置している。この異常データは、機械装置１００８が正常であれば、１％程度の確率でしか生じないとされるデータと判断される。そのため、機械装置１００８に異常が発生している可能性があると判断される。一方、主成分分析を用いずに単純に特徴量に対してしきい値を設けた場合、異常データが正常データの分布に収まってしまい、診断が困難となる場合がある。
　実施の形態５の異常度計算部５４０３は、機械装置１００８の異常の度合いを定量的に示す指標である異常度を計算する。異常度の例として、Ｔ２統計量と呼ばれる指標を用いる。Ｔ２統計量は、主成分分析による写像後のデータが各主成分の軸の中心座標からの距離を軸ごとの分散で標準化した指標である。各主成分の軸の中心座標からの距離は「マハラノビス距離」とも呼ばれる。
　図２０は、実施の形態５における主成分分析の異常度を示す図である。図２０には、図１９に示した正常データ、異常データ、９０％信頼区間、９９％信頼区間を前述の固有値を用いて標準化した場合の図が示されている。また、図２０には、第１主成分の軸と第２主成分の軸とが交差する点と、異常データとの間の距離が異常度として示されている。また、図１９では楕円状となっていた９０％信頼区間及び９９％信頼区間が、図２０では円形になっている。図２０のように、マハラノビス距離の概念を用いて図示すれば、正常時のデータのばらつきを考慮した異常度合いの判断が可能となる。
　実施の形態５では、異常度計算部５４０３として、主成分分析で求めた主成分に写像後のＴ２統計量を用いて異常度を算出する方法を述べたが、この方法に限定されない。Ｔ２統計量の代わりに、Ｑ統計量と呼ばれる指標を用いてもよい。また、主成分分析と同じように正常データを学習してその分布を覚え、新たなデータに対して異常度を算出可能な方法として、教師なし学習と呼ばれる様々な手法が存在する。例えば、１クラスサポートベクターマシン（Ｏｎｅ　ｃｌａｓｓ　Ｓｕｐｐｏｒｔ　Ｖｅｃｔｏｒ　Ｍａｃｈｉｎｅ）、マハラノビス・タグチ法（Ｍａｈａｌａｎｏｂｉｓ　Ｔａｇｕｃｈｉ　Ｍｅｔｈｏｄ）、自己組織化写像（Ｓｅｌｆ－Ｏｒｇａｎｉｚｉｎｇ　Ｍａｐｓ）等の方法を用いてもよい。また、必要に応じて、入力される状態量に正規化処理等の前処理を施してから異常度を計算してもよい。
　実施の形態５の意思決定部５４０４は、前述の方法で計算した異常度を用いて機械装置１００８の状態を診断した結果を決定する。意思決定部５４０４が、特徴量を用いて機械装置１００８の状態を診断する手順について、図２１を参照して説明する。図２１は、実施の形態５における異常度の経年変化の様子の一例を示す図である。図２１には、前述の方法で計算した異常度を定期的にプロットした結果が示されている。定期的とは、例えば一月毎である。図２１の横軸は時間であり、機械装置１００８の運用時間を表している。
　意思決定部５４０４は、異常度に対して、しきい値Ｆｔｈｂを設け、異常度がしきい値Ｆｔｈｂを超えた場合に、機械装置１００８に異常が生じていると診断する。
　図２１において、時刻Ｔｔｂ０は運用を開始した時刻である。時刻Ｔｔｂ０から時刻Ｔｔｂ１より前までの異常度は、多少ばらつきを生じながらも、しきい値Ｆｔｈｂよりも小さい値を保っている。意思決定部５４０４は、異常度がしきい値を下回っている場合には、機械装置１００８が正常である旨の診断を下し、その診断結果を出力する。
　また、図２１の例では、時間の経過と共に異常度の平均的な値が徐々に増加していき、時刻Ｔｔｂ１で異常度がしきい値Ｆｔｈｂを上回っている。この場合、状態診断部５００６は、時刻Ｔｔｂ１で機械装置１００８に異常が生じたという診断を下し、その診断結果を出力する。
　しきい値Ｆｔｈｂの決め方は様々である。例えば正常データの９９％信頼区間に相当する異常度をしきい値Ｆｔｈｂとすることが考えられる。また、過去の別の機械装置に異常が生じた際の異常度をしきい値Ｆｔｈｂとしてもよい。また、機械装置１００８の運用を開始した直後の異常度を参考にして定めてもよい。また、しきい値Ｆｔｈｂは、動的に設定してもよい。例えば複数の機械装置を運用しながら、それらの異常度を個別に計算し、装置ごとの異常度のばらつきを考慮して定期的にしきい値Ｆｔｈｂを決定してもよい。そのように情報処理装置を構成した場合、複数の機械装置のうち一つが故障した際に、故障した装置の異常度がその他の装置の異常度に比べて大きくなるため、異常の検知が可能となる。
　また、本実施の形態では、異常度が一度でもしきい値Ｆｔｈｂを超えた場合に、機械装置１００８が異常であるという診断を下しているが、この手法に限定されない。ノイズ等の影響による誤った診断、即ち誤判定を減らすため、一定の期間の異常度の分布に基づいた検定を行うなどの統計的な手法を採用してもよい。
　以上説明したように、実施の形態５における情報処理装置は、機械装置の異常の度合いを定量的に示す指標である異常度に基づいて機械装置の状態を診断する。これにより、他の実施の形態の効果を得た上で、更に機械装置の状態を高精度に推定できるという効果が得られる。
実施の形態６．
　図２２は、実施の形態６における情報処理装置６０００を含む情報処理システム１００Ｅの構成例を示す図である。図２２に示す情報処理システム１００Ｅでは、図１に示す情報処理システム１００の構成において、状態診断部１００６が状態診断部６００６に置き替えられている。その他の構成は、図１に示す情報処理システム１００と同一又は同等である。同一又は同等の構成部には同一の符号を付すと共に、重複する説明は省略する。なお、本実施の形態では、状態診断部６００６を図１に示す情報処理装置１０００に適用する例で説明するが、図９、図１１又は図１３に示す、情報処理装置２０００，３０００，４０００のうちの何れに適用してもよい。
　図２３は、実施の形態６における状態診断部６００６の構成例を示す図である。状態診断部６００６は、状態量取得部６４０１と、学習部６４０２と、状態推定部６４０３と、意思決定部６４０４とを備える。状態量取得部６４０１は、特徴量を含む状態量を取得する。学習部６４０２は、機械装置１００８の状態と、状態量取得部６４０１が取得する状態量との関係を紐付ける学習用データに基づいて機械装置１００８の状態と当該状態量との関係を学習する。状態推定部６４０３は、学習部６４０２に入力された状態量と学習部６４０２による学習結果とに基づいて機械装置１００８の状態を推定する。意思決定部６４０４は、状態推定部６４０３が推定した結果に基づいて機械装置１００８の状態を診断した診断結果を決定する。
　状態量取得部６４０１は、特徴量と共に、機械装置１００８に付与する設定値を含む情報を状態量として取得する。特徴量の一例として、前述の方法で逐次計算した歪度及び尖度を用いることが考えられる。また、状態量の一つとして用いる機械装置１００８の情報として、前述のモータ速度の設定値を利用することが考えられる。なお、本実施の形態では、状態量取得部６４０１への入力情報が、歪度、尖度及びモータ速度の設定値という３つの場合を例示するが、この例に限定されない。状態量取得部６４０１への入力情報は、４以上であってもよい。
　学習部６４０２及び状態推定部６４０３は、例えば、ニューラルネットワークモデルに従って、いわゆる教師あり学習によって、状態量と機械装置１００８の状態との関係を学習してもよい。ここで、ある入力（状態量）と結果（ラベル）のデータの組である学習用データを大量に学習部に与えることで、それらのデータにある特徴を学習し、入力から結果を推定するモデルを教師あり学習と呼んでいる。ニューラルネットワークは、複数のニューロンからなる入力層、複数のニューロンからなる中間層（隠れ層）、及び複数のニューロンからなる出力層で構成される。中間層は、１層でもよく、２層以上でもよい。
　図２４は、実施の形態６における状態推定部６４０３の構造の例を示す図である。入力は状態量である。説明をわかりやすくするため、図２４のニューラルネットワークは、入力数を３、層数を３としている。複数の入力が、Ｘ１からＸ３で構成される入力層に入力されると、入力値にＷ１１からＷ１６で構成される重みＷ１を乗じた値が、Ｙ１とＹ２とで構成される中間層に入力される。更に、中間層の入力値に、Ｗ２１からＷ２６で構成される重みＷ２を乗じた値が、Ｚ１からＺ３で構成される出力層から出力される。この出力結果は、重みＷ１の値と重みＷ２の値に依存して変化する。
　実施の形態６のニューラルネットワークによる学習の一例では、入力層Ｘ１～Ｘ３に状態量１～３を入力して、出力層Ｚ１～Ｚ３から出力される各状態の確率が学習用データのラベルと一致するように重みＷ１及び重みＷ２を調整する。なお、実施の形態６における学習処理を実行した学習済みの状態推定部６４０３を搭載した情報処理装置６０００を構成してもよい。学習済みの状態推定部６４０３は、学習済のデータ、学習済のプログラム、又はこれらの組み合わせで構成してもよい。学習済みの状態推定部６４０３を用いることにより、他の情報処理装置を用いた学習の結果を利用することができるため、新たに学習を行わずに、診断を実現できる情報処理装置６０００を提供することができる。
　図２４には、歪度が入力層Ｘ１に入力され、尖度が入力層Ｘ２に入力され、モータ速度の設定値が入力層Ｘ３に入力される例が示されている。歪度は入力層Ｘ１に入力される状態量１の例示であり、尖度は入力層Ｘ２に入力される状態量２の例示であり、モータ速度の設定値は入力層Ｘ３に入力される状態量３の例示である。また、出力層Ｚ１からは、状態１の確率として機械装置１００８のボールねじ軸１２２４が故障している確率が出力される。出力層Ｚ２からは、状態２の確率として機械装置１００８のカップリング１２２０が故障している確率が出力される。出力層Ｚ３からは、状態３の確率として機械装置１００８が正常な状態である確率が出力される。出力の計算には、例えば全ての出力の和が１、即ち１００％となるように、例えばソフトマックス関数（ｓｏｆｔｍａｘ　ｆｕｎｃｔｉｏｎ）を用いてもよい。この関数を用いると、推定の結果が理解しやすくなる。
　図２４に示すニューラルネットワークは、学習部６４０２に入力される学習用データに基づいて作成されるデータセットに従って、教師あり学習により、機械装置１００８の状態と状態量１～３との関係を学習する。学習時には、ラベルと呼ばれる機械装置１００８の真の状態を示す情報が正しく得られているものとし、状態１～３のうちの何れかの出力が１（確率が１００％）、それ以外の出力が０（確率が０％）となるように学習するのがよい。
　学習部６４０２及び状態推定部６４０３の構成はニューラルネットワークに限らなくてもよい。例えば、ｋ近傍法、二分木探索、サポートベクターマシン、線形回帰、ロジスティック回帰等、状態量とラベルとの関係を学習し、状態量から状態を推定する方法が数多く知られている。従って、これの手法を適用して、学習部６４０２及び状態推定部６４０３を構成してもよい。
　意思決定部６４０４は、状態推定部６４０３の推定した結果に基づいて機械装置１００８の状態を診断した診断結果を決定する。例えば、機械装置１００８が正常な可能性が最も高い場合、機械装置１００８の診断結果として正常である旨を示す診断結果を出力するのがよい。また、ボールねじ軸１２２４が故障している確率又はカップリング１２２０が故障している確率が最も高い場合、故障箇所を示す診断結果を出力するのがよい。また、ボールねじ軸１２２４が故障している確率とカップリング１２２０が故障している確率との和が、正常な確率よりも大きい場合、何れかの部品が故障していることを示す結果を出力するのがよい。
　以上説明したように、実施の形態６における情報処理装置は、機械装置の状態と状態量との関係を紐付ける学習用データに基づいて機械装置の状態と当該状態量との関係を学習し、その学習結果と学習時に用いた状態量とに基づいて機械装置の状態を推定する。これにより、実施の形態１の効果を得た上で、更に機械装置の状態を高精度に推定できるという効果が得られる。
　以上の実施の形態に示した構成は、一例を示すものであり、別の公知の技術と組み合わせることも可能であるし、実施の形態同士を組み合わせることも可能であるし、要旨を逸脱しない範囲で、構成の一部を省略、変更することも可能である。
　１００，１００Ａ，１００Ｂ，１００Ｃ，１００Ｄ，１００Ｅ　情報処理システム、１０００，２０００，３０００，４０００，５０００，６０００　情報処理装置、１００１，２００１，３００１　センサデータ取得部、１００２　内部変数保持部、１００３，２００３，４００３　内部変数計算部、１００４，２００４，４００４　特徴量計算部、１００５，３００５　初期化処理部、１００６，４００６，５００６，６００６　状態診断部、１００８　機械装置、１００９　モータ、１０１０　センサ、１０９９，２０９９　装置制御部、１２１０　ボールねじ、１２１２　可動部、１２１３　ガイド、１２２０　カップリング、１２２４　ボールねじ軸、１２３０　サーボモータ、１２３１　サーボモータ軸、１２３２　電流センサ、１２３３　エンコーダ、１２４０　ドライバ、１２５０，１２８０　表示器、１２６０　ＰＬＣ、１２７０　ＰＣ、１２９１　プロセッサ、１２９２　メモリ、１２９３　処理回路、４００７　収束度計算部、５４０１，６４０１　状態量取得部、５４０２，６４０２　学習部、５４０３　異常度計算部、６４０３　状態推定部、５４０４，６４０４　意思決定部。

Claims (12)

  1. 　センサによって計測された機械装置の物理量の計測値を取得し、前記計測値のうちの第１時点から第Ｎ時点（Ｎは２以上の整数）までの計測値をセンサデータとして保持するセンサデータ取得部と、
    　前記センサデータに基づいて時系列に順次計算される、前記Ｎよりも少ない個数の内部変数を保持する内部変数保持部と、
    　第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する前記内部変数を、前記第ｊ＋１時点の前記センサデータと第ｊ時点に対応する前記内部変数とに基づいて計算する内部変数計算部と、
    　前記第Ｎ時点の前記内部変数に基づいて前記第１時点から前記第Ｎ時点までの前記センサデータに含まれる統計的な特徴を抽出した特徴量を計算する特徴量計算部と、
    　前記特徴量に基づいて前記機械装置の状態を診断する状態診断部と、
    　を備える情報処理装置。
  2. 　前記特徴量は、前記センサデータの分散、標準偏差、二乗平均平方根、歪度、及び尖度のうちの少なくとも一つである
    　請求項１に記載の情報処理装置。
  3. 　前記機械装置はモータによって駆動され、
    　前記センサデータ取得部は、前記モータの駆動に伴って前記センサで計測される、位置、速度、加速度、動作指令、電流、電圧、トルク、力、圧力、音声、及び光量の計測値のうちの少なくとも一つを前記センサデータとして取得する
    　請求項１又は２に記載の情報処理装置。
  4. 　前記第１時点の前記内部変数を予め設定した最大値と予め設定した最小値との間の値に決定する初期化処理を実行する初期化処理部を備える
    　請求項１から３の何れか１項に記載の情報処理装置。
  5. 　前記センサデータ取得部は、前記センサデータを取得する時間の間隔である取得周期、前記機械装置の動作の状況を表す動作信号、及び前記センサデータにおけるデータ値のうちの少なくとも一つに基づいて前記第１時点から前記第ｊ時点までの間の２点以上の時点を決定し、決定した前記２点以上の時点に基づいて内部変数を更新するための計算許可フラグを生成し、
    　前記内部変数計算部は、前記計算許可フラグの示す各第ｕ時点（ｕは１以上、且つ前記ｊ以下の整数）について、第ｕ＋１時点の前記内部変数を、前記第ｕ時点の前記内部変数の値に決定する
    　請求項４に記載の情報処理装置。
  6. 　前記センサデータ取得部は、前記センサデータを取得する時間の間隔である取得周期、前記機械装置の動作の状況を表す動作信号、及び前記センサデータのうちの少なくとも一つに基づいて前記第１時点から第Ｎ時点までの間の１点以上の時点を決定し、決定した前記１点以上の時点に基づいて内部変数を初期化するための初期化トリガを生成し、
    　前記初期化処理部は、前記初期化トリガに基づいて前記内部変数を初期化する処理を実行する
    　請求項４に記載の情報処理装置。
  7. 　前記内部変数計算部は、時刻が古い前記センサデータが前記特徴量に与える影響よりも時刻が新しい前記センサデータのほうが前記特徴量に与える影響を大きくするように重みをつけるための、０より大きく１より小さい忘却係数に基づいて前記内部変数を計算する
    　請求項１から６の何れか１項に記載の情報処理装置。
  8. 　前記内部変数に基づいて、前記特徴量の計算の収束の度合いを定量的に示す指標である収束度を計算する収束度計算部を備える
    　請求項７に記載の情報処理装置。
  9. 　前記状態診断部は、
    　前記特徴量を含む状態量を取得する状態量取得部と、
    　前記機械装置が正常な状態の状態量に基づいて前記機械装置の状態と前記特徴量との関係を学習する学習部と、
    　前記学習部が学習した結果に基づいて前記機械装置の異常の度合いを定量的に示す指標である異常度を計算する異常度計算部と、
    　前記異常度に基づいて前記機械装置の状態を診断した診断結果を決定する意思決定部と、
    　を備える請求項１から８の何れか１項に記載の情報処理装置。
  10. 　前記状態診断部は、
    　前記特徴量を含む状態量を取得する状態量取得部と、
    　前記機械装置の状態と前記状態量との関係を紐付ける学習用データに基づいて前記機械装置の状態と前記状態量との関係を学習する学習部と、
    　前記学習部に入力された状態量と前記学習部による学習結果とに基づいて前記機械装置の状態を推定する状態推定部と、
    　前記状態推定部が推定した結果に基づいて前記診断の結果を決定する意思決定部と、
    　を備える請求項１から９の何れか１項に記載の情報処理装置。
  11. 　機械装置の物理量を計測するセンサから、第１時点から第Ｎ時点（Ｎは２以上の整数）までの各時点における前記物理量の計測値をセンサデータとして取得する第１ステップと、
    　前記センサデータに基づいて時系列に順次計算される、前記Ｎよりも少ない個数の内部変数を保持する第２ステップと、
    　第ｊ＋１時点（ｊは１からＮ－１までの整数）に対応する前記内部変数を、前記第ｊ＋１時点の前記センサデータと第ｊ時点に対応する前記内部変数とに基づいて計算する第３ステップと、
    　前記第Ｎ時点の前記内部変数に基づいて前記第１時点から前記第Ｎ時点までの前記センサデータに含まれる統計的な特徴を抽出した特徴量を計算する第４ステップと、
    　前記特徴量に基づいて前記機械装置の状態を診断する第５ステップと、
    　を含む情報処理方法。
  12. 　前記第１時点の前記内部変数を予め設定した最大値と最小値との間の値に決定する初期化処理を実行するステップを更に有する
    　請求項１１に記載の情報処理方法。
PCT/JP2020/047394 2020-12-18 2020-12-18 情報処理装置及び情報処理方法 Ceased WO2022130611A1 (ja)

Priority Applications (7)

Application Number Priority Date Filing Date Title
US18/031,172 US12498709B2 (en) 2020-12-18 2020-12-18 Information processing apparatus and information processing method
KR1020237018898A KR102841817B1 (ko) 2020-12-18 2020-12-18 정보 처리 장치 및 정보 처리 방법
JP2021526355A JP6935046B1 (ja) 2020-12-18 2020-12-18 情報処理装置及び情報処理方法
DE112020007851.5T DE112020007851T5 (de) 2020-12-18 2020-12-18 Informationsverarbeitungsvorrichtung und informationsverarbeitungsverfahren
PCT/JP2020/047394 WO2022130611A1 (ja) 2020-12-18 2020-12-18 情報処理装置及び情報処理方法
CN202080107328.7A CN116569120A (zh) 2020-12-18 2020-12-18 信息处理装置及信息处理方法
TW110146274A TWI814170B (zh) 2020-12-18 2021-12-10 資訊處理裝置以及資訊處理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2020/047394 WO2022130611A1 (ja) 2020-12-18 2020-12-18 情報処理装置及び情報処理方法

Publications (1)

Publication Number Publication Date
WO2022130611A1 true WO2022130611A1 (ja) 2022-06-23

Family

ID=77657855

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/047394 Ceased WO2022130611A1 (ja) 2020-12-18 2020-12-18 情報処理装置及び情報処理方法

Country Status (7)

Country Link
US (1) US12498709B2 (ja)
JP (1) JP6935046B1 (ja)
KR (1) KR102841817B1 (ja)
CN (1) CN116569120A (ja)
DE (1) DE112020007851T5 (ja)
TW (1) TWI814170B (ja)
WO (1) WO2022130611A1 (ja)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2023051402A (ja) * 2021-09-30 2023-04-11 オムロン株式会社 制御システム、情報処理方法および情報処理装置
JP2023065122A (ja) * 2021-10-27 2023-05-12 株式会社ジェイテクト 機械設備及びボールねじ機構の診断方法
JP7638454B2 (ja) * 2022-09-08 2025-03-03 三菱電機株式会社 異常要因推定装置、学習装置、精密診断システム、および、異常要因推定方法
CN115859201B (zh) * 2022-11-22 2023-06-30 淮阴工学院 一种化工过程故障诊断方法及系统
JPWO2024134843A1 (ja) * 2022-12-22 2024-06-27
WO2024142315A1 (ja) * 2022-12-27 2024-07-04 三菱電機株式会社 機器診断装置、プログラム、機器診断システム及び機器診断方法
US12594938B2 (en) * 2023-06-21 2026-04-07 Dana Heavy Vehicle Systems Group, Llc Optimized electric machine stop position for loss reduction and vehicle launch
KR102834176B1 (ko) * 2023-11-17 2025-07-15 한국공학대학교산학협력단 감속기의 백래시 평가 모델을 생성하는 방법, 감속기의 백래시를 보상하는 방법 및 장치

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH05187965A (ja) * 1991-08-12 1993-07-27 Kuroda Precision Ind Ltd ボールねじの寿命監視装置
JP2002341909A (ja) * 2001-05-18 2002-11-29 Sofutorokkusu:Kk 工作機器の監視方法
CN102033200A (zh) * 2009-09-29 2011-04-27 上海宝钢工业检测公司 基于统计模型的交流电机在线监测和诊断方法
WO2013105164A1 (ja) * 2012-01-13 2013-07-18 日本電気株式会社 異常信号判定装置、異常信号判定方法、および異常信号判定プログラム
JP2017162417A (ja) * 2016-03-11 2017-09-14 株式会社日立ハイテクノロジーズ 異常診断装置および方法、並びに、異常診断システム
JP2020035372A (ja) * 2018-08-31 2020-03-05 オムロン株式会社 情報処理装置及び情報処理方法

Family Cites Families (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2603134B2 (ja) * 1989-06-06 1997-04-23 三菱電機株式会社 移動平均処理装置
DE10034224B4 (de) * 2000-07-07 2016-11-17 Continental Teves Ag & Co. Ohg Verfahren zur Bildung eines Mittelwertes bei der Erkennung eines Druckverlustes in Kraftfahrzeugreifen
JP2003329493A (ja) * 2002-05-14 2003-11-19 Japan Research Institute Ltd 状態監視方法、状態監視システム、状態監視装置、コンピュータプログラム、及び記録媒体
WO2004063939A2 (de) 2003-01-10 2004-07-29 Robert Bosch Gmbh Verfahren zum berechnen eines mittelwertes von messwerten
US7075327B2 (en) * 2003-06-18 2006-07-11 Eaton Corporation System and method for proactive motor wellness diagnosis
JP2008261755A (ja) 2007-04-12 2008-10-30 Canon Inc 情報処理装置、情報処理方法
WO2012009804A1 (en) * 2010-07-23 2012-01-26 Corporation De L'ecole Polytechnique Tool and method for fault detection of devices by condition based maintenance
JP2013191113A (ja) 2012-03-15 2013-09-26 Sony Corp 表示制御装置、表示制御方法およびプログラム
JP5485441B2 (ja) * 2013-03-19 2014-05-07 株式会社日立製作所 異常診断装置および産業機械
JP5771317B1 (ja) * 2014-08-26 2015-08-26 株式会社日立パワーソリューションズ 異常診断装置及び異常診断方法
US11327475B2 (en) * 2016-05-09 2022-05-10 Strong Force Iot Portfolio 2016, Llc Methods and systems for intelligent collection and analysis of vehicle data
WO2018003457A1 (ja) * 2016-06-30 2018-01-04 パナソニックＩｐマネジメント株式会社 情報処理装置、時系列データの情報処理方法、及びプログラム
JP6840953B2 (ja) * 2016-08-09 2021-03-10 株式会社リコー 診断装置、学習装置および診断システム
JP6723669B2 (ja) * 2016-09-27 2020-07-15 東京エレクトロン株式会社 異常検知プログラム、異常検知方法および異常検知装置
JP6842299B2 (ja) * 2016-12-28 2021-03-17 三菱パワー株式会社 診断装置、診断方法及びプログラム
JP7085370B2 (ja) * 2017-03-16 2022-06-16 株式会社リコー 診断装置、診断システム、診断方法およびプログラム
DE102017217804A1 (de) * 2017-10-06 2019-04-11 Bayerische Motoren Werke Aktiengesellschaft Vorrichtung und verfahren zur ermittlung einer tachometerkennlinie eines fahrzeugs, system zur regelung der geschwindigkeit eines fahrzeugs sowie fahrzeug
JP7283485B2 (ja) * 2018-12-28 2023-05-30 日本電気株式会社 推定装置、推定方法、及びプログラム

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH05187965A (ja) * 1991-08-12 1993-07-27 Kuroda Precision Ind Ltd ボールねじの寿命監視装置
JP2002341909A (ja) * 2001-05-18 2002-11-29 Sofutorokkusu:Kk 工作機器の監視方法
CN102033200A (zh) * 2009-09-29 2011-04-27 上海宝钢工业检测公司 基于统计模型的交流电机在线监测和诊断方法
WO2013105164A1 (ja) * 2012-01-13 2013-07-18 日本電気株式会社 異常信号判定装置、異常信号判定方法、および異常信号判定プログラム
JP2017162417A (ja) * 2016-03-11 2017-09-14 株式会社日立ハイテクノロジーズ 異常診断装置および方法、並びに、異常診断システム
JP2020035372A (ja) * 2018-08-31 2020-03-05 オムロン株式会社 情報処理装置及び情報処理方法

Also Published As

Publication number Publication date
CN116569120A (zh) 2023-08-08
TW202225886A (zh) 2022-07-01
US12498709B2 (en) 2025-12-16
JPWO2022130611A1 (ja) 2022-06-23
JP6935046B1 (ja) 2021-09-15
US20230376023A1 (en) 2023-11-23
TWI814170B (zh) 2023-09-01
KR102841817B1 (ko) 2025-08-01
DE112020007851T5 (de) 2023-09-28
KR20230104236A (ko) 2023-07-07

Similar Documents

Publication Publication Date Title
JP6935046B1 (ja) 情報処理装置及び情報処理方法
JP6810097B2 (ja) 異常検出器
CN110187694B (zh) 故障预测装置以及机器学习装置
CN106409120B (zh) 机械学习方法及机械学习装置、以及故障预知装置及系统
US11235461B2 (en) Controller and machine learning device
JP6451662B2 (ja) 異常判定装置、異常判定プログラム、異常判定システム、及びモータ制御装置
US10684608B2 (en) Abnormality detection apparatus and machine learning device
CN113219934B (zh) 自动化系统中的部件的参数设定
US12103169B2 (en) Abnormality diagnosis device and abnormality diagnosis method
CN117222878A (zh) 机器人故障预兆检测装置以及机器人故障预兆检测方法
US10747197B2 (en) Abnormally factor identification apparatus
US10802476B2 (en) Numerical controller with learned pressure estimation
WO2020213062A1 (ja) モータ制御装置
JP2021096639A (ja) 制御方法、制御装置、機械設備、制御プログラム、記録媒体
CN116635802B (zh) 数控装置
US10962957B2 (en) Collision position estimation device and machine learning device
JP2019024305A (ja) モータ制御システム
WO2021019760A1 (ja) 異常診断方法、異常診断装置および異常診断プログラム
JPWO2020152741A1 (ja) 異常要因推定装置、異常要因推定方法、及びプログラム
JP7011106B1 (ja) 状態判定装置及び状態判定方法
US20230191513A1 (en) Tool diagnostic device
CN114693143A (zh) 一种数控机床的健康状态评价方法、系统、设备和介质
JP2022136350A (ja) 時系列データ解析に基づいたシステム制御
US12596354B2 (en) Abnormality information estimation system, operation analysis system, motor control device, abnormality information estimation method, and program
WO2024176390A1 (ja) データ分割装置、及びコンピュータ読み取り可能な記録媒体

Legal Events

Date Code Title Description
ENP Entry into the national phase

Ref document number: 2021526355

Country of ref document: JP

Kind code of ref document: A

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

Ref document number: 20965998

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 202327029905

Country of ref document: IN

WWE Wipo information: entry into national phase

Ref document number: 202080107328.7

Country of ref document: CN

ENP Entry into the national phase

Ref document number: 20237018898

Country of ref document: KR

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 112020007851

Country of ref document: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20965998

Country of ref document: EP

Kind code of ref document: A1

WWG Wipo information: grant in national office

Ref document number: 202327029905

Country of ref document: IN

WWG Wipo information: grant in national office

Ref document number: 1020237018898

Country of ref document: KR

WWG Wipo information: grant in national office

Ref document number: 18031172

Country of ref document: US