WO2019163701A1 - システム同定装置、システム同定方法及び記録媒体 - Google Patents

システム同定装置、システム同定方法及び記録媒体 Download PDF

Info

Publication number
WO2019163701A1
WO2019163701A1 PCT/JP2019/005805 JP2019005805W WO2019163701A1 WO 2019163701 A1 WO2019163701 A1 WO 2019163701A1 JP 2019005805 W JP2019005805 W JP 2019005805W WO 2019163701 A1 WO2019163701 A1 WO 2019163701A1
Authority
WO
WIPO (PCT)
Prior art keywords
system identification
response function
unit
analysis
analysis target
Prior art date
Application number
PCT/JP2019/005805
Other languages
English (en)
French (fr)
Inventor
宗一朗 高田
裕文 井上
茂樹 篠田
菊池　克
Original Assignee
日本電気株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 日本電気株式会社 filed Critical 日本電気株式会社
Priority to EP19756575.7A priority Critical patent/EP3757704A4/en
Priority to JP2020501746A priority patent/JP6981526B2/ja
Priority to US16/970,754 priority patent/US20210010980A1/en
Publication of WO2019163701A1 publication Critical patent/WO2019163701A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4472—Mathematical theories or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M7/00—Vibration-testing of structures; Shock-testing of structures
    • G01M7/02—Vibration-testing by means of a shake table
    • G01M7/025—Measuring arrangements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H13/00—Measuring resonant frequency
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M7/00—Vibration-testing of structures; Shock-testing of structures
    • G01M7/08—Shock-testing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G01N29/12—Analysing solids by measuring frequency or resonance of acoustic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4409—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison
    • G01N29/4418—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison with a model, e.g. best-fit, regression analysis
    • 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
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/26—Scanned objects
    • G01N2291/262—Linear objects
    • G01N2291/2626—Wires, bars, rods

Definitions

  • the present invention relates to a system identification device, a system identification method, and a recording medium.
  • the system to be modeled is a system having proximity eigenvalues such as a beat phenomenon or resonance
  • system identification may be difficult with the techniques described in Patent Documents 1 to 4 and Non-Patent Document 1.
  • the system having the proximity eigenvalue includes a system having a multiple root in which the eigenvalue is degenerated.
  • An object of this invention is to provide the system identification apparatus etc. which can perform the system identification of the system which has a proximity
  • the system identification device calculates a self-frequency response function based on an input signal and an output signal measured at a position where the analysis target is excited, and the calculated self-frequency response
  • An analysis unit that performs system identification of the analysis target using an impulse response function obtained from the function and an impulse response function of a virtual two-degree-of-freedom model in which the analysis target is modeled.
  • the system identification method includes an excitation step for exciting an analysis target, and a measurement step for measuring an input signal and an output signal at a position where the analysis target is excited in the excitation step.
  • a self-frequency response function is calculated based on the input signal and the output signal measured in the measurement step, and the impulse response function obtained from the calculated self-frequency response function and the analysis target are modeled.
  • the recording medium calculates a self-frequency response function based on an input signal and an output signal measured at a position where the analysis target is excited by the computer, and the calculated self
  • a program for executing an analysis step of performing system identification of the analysis target using an impulse response function obtained from a frequency response function and an impulse response function of a virtual two-degree-of-freedom model in which the analysis target is modeled is recorded.
  • system identification of a system having proximity eigenvalues can be performed.
  • system identification problem The mathematical modeling method of physical systems based on observation data is called “system identification problem”. This problem is broadly classified into (1) a case where the input and output signals of the system are known and (2) a case where the input is unknown. Furthermore, a technique using a time domain signal or a frequency domain signal is also known.
  • AR model Autoregressive model: autoregressive model
  • MA model Moving average model
  • ARMA model Autoregressive moving average model
  • ARX model A polynomial model such as Auto-Regressive (eXogeneous model) is used.
  • the polynomial model is difficult to apply to a system having adjacent eigenvalues in which the frequency domain is flattened and different eigenvalues are close to each other.
  • frequency domain information when frequency domain information is used, the peak positions of both eigenvalues are often unclear, and the curve fitting method cannot be applied. In particular, the peak position cannot be visually recognized depending on the number of samples in the frequency domain by Fourier transform.
  • the system identification problem of the system having the proximity eigenvalue is a problem that has not been sufficiently solved.
  • the system identification apparatus or the like according to the present embodiment solves such a problem and performs system identification of a system having a proximity eigenvalue (including a system having a multiple root where the eigenvalues overlap).
  • FIG. 1 is a block diagram showing a configuration of a system identification device 1 according to the first embodiment.
  • the system identification device 1 includes an installation positioning unit 101, an excitation unit 102, a measurement unit 103, a signal collection unit 104, and an analysis unit 105.
  • the target physical system 106 is an identification target by the system identification device 1.
  • the installation positioning unit 101 installs the excitation unit 102 and the measurement unit 103 in the target physical system 106.
  • the vibration unit 102 excites the target physical system 106 via the installation positioning unit 101.
  • the measurement unit 103 detects an input signal and an output signal to the target physical system 106 when the excitation unit 102 excites the target physical system 106 via the installation positioning unit 101.
  • the signal collection unit 104 converts the input signal and output signal detected by the measurement unit 103 into data.
  • the analysis unit 105 analyzes the data obtained by the signal collection unit 104 and performs system identification of the target physical system 106.
  • FIG. 2 is a flowchart showing processing of the system identification device 1.
  • the measurer installs the excitation unit 102 and the measurement unit 103 via the installation positioning unit 101 in the target physical system 106 (step S110).
  • the input position is matched with the output position.
  • the input position is a position where the excitation unit 102 is installed, that is, a cause of vibration is input.
  • the output position is a position where the measurement unit 103 is set, that is, a position where vibration of the target physical system 106 is measured.
  • the target physical system 106 is excited by the vibration unit 102 via the installation positioning unit 101.
  • the measurement unit 103 detects a vibration input signal to the target physical system 106 and a vibration output signal from the target physical system 106.
  • the signal collection unit 104 converts the input signal and output signal detected by the measurement unit 103 into data, and outputs the data to the analysis unit 105.
  • the analysis unit 105 analyzes the obtained data. Specifically, the analysis unit 105 first performs fast Fourier transform (FFT) on the input signal and the output signal, respectively.
  • FFT fast Fourier transform
  • the analysis unit 105 obtains a self-frequency response function by dividing the output signal by the input signal in the frequency domain (step S120).
  • the analysis unit 105 performs zooming only in the frequency band where the target eigenvalue exists in the self-frequency response function (step S130).
  • the analyzing unit 105 obtains an impulse response function of the self-frequency response function by performing an inverse Fourier transform on the zoomed self-frequency response function (Step S140).
  • the analysis unit 105 receives an initial value and a step width used in the next step S160 (step S150).
  • the analysis unit 105 applies the multivariable Newton method using the impulse response function of the virtual two-degree-of-freedom system described later to the impulse response function obtained in step S140 (step S160).
  • the analysis unit 105 uses the initial value and step width input in step S150 when the multivariable Newton method is executed.
  • step S170 determines that the solution does not converge (step S170: NO)
  • the analysis unit 105 returns to step S150, receives the input of the initial value and step width to be newly used, and performs the process of step S160.
  • step S170: YES the mass, stiffness constant, and attenuation of the system of the target physical system 106 based on the impulse response function of the virtual two-degree-of-freedom system at the time of convergence.
  • a coefficient is obtained and system identification is performed (step S180). According to this embodiment, system identification of a system having a proximity eigenvalue is possible.
  • FIG. 3 is a block diagram showing the configuration of the system identification device 3 according to the second embodiment.
  • the system identification device 3 includes an installation positioning unit 301, an excitation unit 302, a measurement unit 303, a signal collection unit 304, an analysis unit 305, a storage unit 306, and an initial value setting unit 307.
  • the target physical system 308 is a pipeline and is an identification target by the system identification device 3.
  • the installation positioning unit 301, the excitation unit 302, the measurement unit 303, and the signal collection unit 304 have the same functions as the installation positioning unit 101, the excitation unit 102, the measurement unit 103, and the signal collection unit 104 of the first embodiment, respectively.
  • the storage unit 306 stores pipeline management ledger data.
  • the pipe management ledger data includes information on the diameter, material type, and pipe thickness, which are physical data of the target physical system 308.
  • the initial value setting unit 307 calculates initial values of the parameters of the multivariable Newton method used in the analysis unit 305 based on the diameter, material type, and tube thickness data read from the storage unit 306.
  • the analysis unit 305 has the same function as the analysis unit 105 of the first embodiment, except that the initial value calculated by the initial value setting unit 307 is used in the multivariable Newton method.
  • FIG. 4 is a flowchart showing the processing of the system identification device 3.
  • the measurer installs the excitation unit 302 and the measurement unit 303 via the installation positioning unit 301 in the target physical system 308 (step S310).
  • the target physical system 308 is excited by the vibration unit 302 via the installation positioning unit 301.
  • the measurement unit 303 detects an input signal to the target physical system 308 and an output signal from the target physical system 308 at the excitation position.
  • the signal collection unit 304 converts the input signal and output signal detected by the measurement unit 303 into data.
  • the analysis unit 305 obtains a self-frequency response function using the input signal and the output signal (step S320).
  • the analysis unit 305 performs zooming only in the frequency band where the target eigenvalue exists in the self-frequency response function (step S330).
  • the analysis unit 305 performs inverse Fourier transform on the zooming portion to obtain an impulse response function of the self-frequency response function (step S340).
  • the initial value setting unit 307 reads data on the diameter, material type, and tube thickness of the target physical system 308 from the storage unit 306 (step S350).
  • the initial value setting unit 307 calculates initial values of the parameters of the multivariable Newton method based on the read data (step S360).
  • the analysis unit 305 receives an input of the step width (step S370).
  • the analysis unit 305 applies the multivariable Newton method using the virtual two-degree-of-freedom impulse response function described later to the impulse response function obtained in step S340 (step S380). In the multivariable Newton method, the initial value calculated in step S360 and the step width input in step S370 are used.
  • step S390 NO
  • the analysis unit 305 determines that the solution does not converge (step S390: NO)
  • the analysis unit 305 returns to step S370, receives the input of the step width to be newly used, and performs the process of step S380.
  • step S390: YES the mass, stiffness constant, and attenuation of the target physical system 308 based on the impulse response function of the virtual two-degree-of-freedom system at the time of convergence.
  • a coefficient is obtained and system identification is performed (step S400). According to this embodiment, the system identification of a pipeline is possible.
  • FIG. 5 is a block diagram showing the configuration of the system identification device 5 of the present embodiment.
  • the system identification device 5 includes a fire hydrant coupler 501, a hammer 502, a sensor 503, a data logger 504, and an identification processing unit 505.
  • the pipe 506 is a water pipe that is a system identification target by the system identification device 5.
  • the hammer 502 include an impulse hammer with a built-in force sensor, a commercially available hammer with an acceleration pickup, and an electromagnetic vibrator.
  • Examples of the sensor 503 include an acceleration pickup, a laser Doppler type speedometer, a laser displacement meter, and a contact displacement meter.
  • the identification processing unit 505 is realized by, for example, a processor, a memory, and an HDD (Hard disk drive).
  • the processor operates as the identification processing unit 505 by reading an identification processing program for causing the computer to execute the processing from the HDD and executing it.
  • FIG. 6 is a diagram showing the state of the fire hydrant coupler 501 installed in the pipe 506.
  • the measurer installs the fire hydrant coupler 501 and the sensor 503 in the pipeline 506 (step S110 in FIG. 2).
  • the measurer strikes the fire hydrant coupler 501 shown in FIG. 6 with the hammer 502 and excites the pipe 506.
  • the sensor 503 detects an output signal after excitation at the same measurement position (Measurements point for vibration response) as the striking position (Tapping Point).
  • the data logger 504 collects the input signal of the hammer 502 and the output signal of the sensor 503.
  • the data logger 504 converts the collected input signal and output signal into data and outputs the data to the identification processing unit 505.
  • the identification processing unit 505 performs a fast Fourier transform (FFT) process on each of the input signal and the output signal.
  • the spectrum of the input signal obtained by FFT is represented as X ( ⁇ ), and the spectrum of the output signal is represented as Y ( ⁇ ).
  • ⁇ is a frequency.
  • L ( ⁇ ) (Y ( ⁇ ) ⁇ X * ( ⁇ )) / (X ( ⁇ ) ⁇ X * ( ⁇ )) is estimated as H 1
  • L ( ⁇ ) (Y ( ⁇ ) ⁇ Y * ( ⁇ )) / (X ( ⁇ ) ⁇ Y * ( ⁇ )) is called H 2 estimation, either of which may be used.
  • X * ( ⁇ ) is the complex conjugate of X ( ⁇ )
  • Y * ( ⁇ ) is the complex conjugate of Y ( ⁇ ).
  • the identification processing unit 505 performs zooming only on the frequency band in which the close eigenvalue of interest exists (step S130 in FIG. 2). A peak appears in the frequency band where the close proximity eigenvalue of interest exists. Therefore, the identification processing unit 505 identifies the frequency band in which the noted eigenvalue is present by detecting a peak in the self-frequency response function L ( ⁇ ). Alternatively, the identification processing unit 505 may display the self-frequency response function L ( ⁇ ) on a display device included in the system identification device 5, and a user who confirms the display may input a frequency band in which a peak appears.
  • the identification processing unit 505 extracts the self-frequency response function L ( ⁇ ) of the specified frequency band, and performs zooming in which a value smaller than the threshold is replaced with 0.
  • the identification processing unit 505 obtains an impulse response function g e (t) by performing inverse FFT on the zoomed result. Note that t represents time (step S140 in FIG. 2).
  • a virtual two-degree-of-freedom model is virtualized as a system identification model.
  • the virtual two-degree-of-freedom model is also called a symmetric two-degree-of-freedom spring mass system.
  • FIG. 7 is a diagram illustrating a virtual two degree of freedom model.
  • the virtual two-degree-of-freedom model is a system in which a one-degree-of-freedom spring mass system having an equal mass, a spring constant, and a damping coefficient is connected by a spring and a dashpot.
  • M is a mass
  • K is a spring constant
  • C is a damping coefficient
  • F is an external force vector
  • x 1 and x 2 are displacement vectors.
  • delta K changes in spring constant
  • delta C shows the variation of the attenuation coefficient.
  • the virtual two-degree-of-freedom model is a system in which the mass matrix, the stiffness matrix, and the damping matrix representing the equation of motion are symmetric matrices, and has the property that the eigenvectors are symmetric, which is suitable for system identification of systems with close eigenvalues. It becomes.
  • ⁇ (t) is a Dirac delta function.
  • the first term (1 / M) ⁇ ⁇ (t) is negligibly smaller than the other terms.
  • each parameter ⁇ d1 , ⁇ d2 , ⁇ , ⁇ , ⁇ is defined as equation (3).
  • An update equation for parameter estimation is obtained by the multivariable Newton method using the square sum J of the difference between the impulse response function g e (t) of the experimental value and the equation (2) as an objective function.
  • the objective function J is shown in Formula (4), and the update formula is shown in Formula (5).
  • i represents the time sampling of the number, is t i went to i-th sampling.
  • ⁇ is a parameter for step adjustment, and when the parameter estimation algorithm diverges, adjusting between 0.001 and 0.1 improves the convergence.
  • ⁇ is a parameter for step adjustment, and when the parameter estimation algorithm diverges, adjusting between 0.001 and 0.1 improves the convergence.
  • g ⁇ 11 is obtained by calculating Expression (2) using the current parameter value.
  • initial values of parameters ( ⁇ d1 , ⁇ d2 , ⁇ , ⁇ , ⁇ ) and a value of ⁇ that is a step width are input.
  • the identification processing unit 505 may receive an input of information used for calculating an initial value and calculate an initial value based on the input information.
  • the identification processing unit 505 is based on the impulse response function g e (t i ) obtained in step S140 and g 11 (t i ) calculated by the equation (2) using the current value of each parameter. Then, the sum of squares J is calculated by equation (4). The identification processing unit 505 updates the value of each parameter according to the equation (5) so that the sum of squares J is equal to or less than the threshold (step S160 in FIG. 2).
  • the identification processing unit 505 repeatedly updates the parameter, and determines whether or not the value has converged due to the value of J or a change in the value. If the identification processing unit 505 determines not to converge (step S170 in FIG. 2: NO), the identification processing unit 505 receives a new initial value and ⁇ (step S150). If the identification processing unit 505 determines that it has converged, the parameter value at that time is used to calculate the mass M, the spring constant K, and the damping coefficient C based on the relationship of Expression (5) (step S180).
  • FIG. 8 is a diagram showing a system identification result under the above conditions.
  • the horizontal axis represents frequency
  • the vertical axis represents acceleration.
  • the figure shows a self-frequency response function (Identified) calculated using the identified parameters and a true self-frequency response function (Experiment). According to the figure, a good match between the true value and the identification result can be confirmed.
  • the initial values of the parameters ⁇ d1 , ⁇ d2 , ⁇ , ⁇ , ⁇ are calculated as follows.
  • the mass M and the spring constant K, which are the main parameters that determine these initial values, are calculated using the following equation (6).
  • R is a radius
  • A is a cross-sectional area.
  • the radius R and the tube thickness h are obtained from the diameter and the tube thickness read from the pipeline management ledger data.
  • the tube length L may be read from the pipeline management ledger data or may be input by the user.
  • E is the modulus of elasticity of the tube
  • ⁇ is the density of the tube.
  • the elastic coefficient E and the pipe density ⁇ may be values corresponding to the material type read from the pipe management ledger data, or may be predetermined values.
  • Delta K is 1/100 of the spring constant K
  • delta C is desirably about 1/3 of the attenuation coefficient C.
  • Attenuation coefficient used for calculating the initial value C, the change delta C of the damping coefficient, diameter read from line management register data, pipe thickness, may be a value corresponding to one or more of grades, a predetermined value It is good.
  • the initial value setting unit 307 uses these values to calculate initial values of the parameters ⁇ d1 , ⁇ d2 , ⁇ , ⁇ , and ⁇ using Expression (3). The above is a specific expression of the initial value setting unit.
  • the in-service water pipe was installed in the test pipeline, and a system identification experiment was conducted in a flowing water environment.
  • the test pipe line used was an ordinary cast iron pipe having a diameter of 100 mm and a wall thickness of 10 mm.
  • a fire hydrant was installed at the top of the pipeline, and an acceleration sensor was installed on the coupler.
  • FIG. 9 is a diagram illustrating the results of a system identification experiment.
  • the horizontal axis represents frequency
  • the vertical axis represents acceleration.
  • the circle plot in the figure represents the experimental value (Experiment) of the self-frequency response function
  • the solid line represents the identification result (Identified).
  • M 14.8446 kg
  • K 3.4346 ⁇ 10 9 N / m
  • C 903.5491 Ns / m
  • ⁇ K 1.2877 ⁇ 10 8 N / m
  • ⁇ C 903.5491 Ns / m.
  • FIG. 10 is a diagram illustrating a comparison result between the identification method using the AR model and the identification method according to the present embodiment.
  • the horizontal axis represents frequency
  • the vertical axis represents acceleration.
  • the circle plot in the figure represents the experimental value (Experiment) of the self-frequency response function
  • the symbol L1 represents the identification result (Identified) of the present embodiment.
  • symbol L2 shows the identification result (AR) by AR method.
  • the calculation condition is that the AR model order is 100th, and the impulse response function of this experiment is used as an identification input. This is exactly the same as the evaluation signal of the system identification method of this embodiment.
  • the identification result by the AR method confirms a large deviation from the experimental value and confirms that identification is difficult. This indicates that the impulse response function is accompanied by a beat waveform, and the polynomial model used in the time domain identification method has a limit in describing the characteristics of this waveform.
  • FIG. 11 is a block diagram showing the minimum configuration of the system identification apparatus according to the embodiment of the present invention.
  • the system identification apparatus 1a with the minimum configuration shown in the figure may include at least the analysis unit 105 of the first embodiment described above.
  • the analysis unit 105 calculates a self-frequency response function based on the input signal and the output signal measured at the position where the analysis target is excited. Then, the analysis unit 105 performs system identification of the analysis target using the impulse response function obtained from the calculated self-frequency response function and the impulse response function of the virtual two-degree-of-freedom model in which the analysis target is modeled. .
  • the system identification apparatuses 1, 1a, 3, and 5 in the above-described embodiment include a CPU (Central Processing Unit) connected via a bus, a memory, an auxiliary storage device, and the like, and the above-described embodiment by executing a system identification program.
  • the functions of the system identification devices 1, 1a, 3, and 5 in FIG. Some functions of the system identification devices 1, 1a, 3, and 5 are realized by using hardware such as ASIC (Application Specific Integrated Circuit), PLD (Programmable Logic Device), and FPGA (Field Programmable Gate Array). May be.
  • the system identification program may be recorded on a computer-readable recording medium.
  • the computer-readable recording medium is a portable medium such as a flexible disk, a magneto-optical disk, a ROM (Read Only Memory), a CD-ROM (Compact Disc Read Only memory), or a hard disk built in a computer system.
  • the system identification program may be transmitted via a telecommunication line.
  • a self-frequency response function is calculated based on an input signal and an output signal measured at a position where the analysis target is excited, an impulse response function obtained from the calculated self-frequency response function, and the analysis target
  • a system identification apparatus comprising: an analysis unit that performs system identification of the analysis target using an impulse response function of a virtual two-degree-of-freedom model in which is modeled.
  • the analysis unit estimates the impulse response function of the virtual two-degree-of-freedom model by a multivariable Newton method, and performs system identification of the analysis target based on the impulse response function obtained by the estimation.
  • the system identification apparatus according to 1.
  • Additional remark 3 The system identification apparatus of Additional remark 2 further provided with the initial value setting part which calculates the initial value used for the multivariable Newton method based on the physical data of the said analysis object.
  • Additional remark 4 The system identification apparatus of Additional remark 1 whose said analysis object is a pipe line.
  • the excitation part which excites the said analysis object The measurement part which measures the input signal and output signal in the position which excited the said analysis object by the said excitation part, The said excitation part is said analysis object
  • the said vibration part is a system identification apparatus of Additional remark 5 which is an impulse hammer or electromagnetic vibrator which incorporated the force sensor.
  • Excitation Step for Exciting Analytical Object Measuring Step for Measuring Input Signal and Output Signal at the Position at which the Analyzing Object is Excited in the Exciting Step, and the Input Signal Measured in the Measuring Step
  • an analysis step for performing system identification of the analysis object using the system identification method

Abstract

システム同定装置１の分析部１０５は、加振部１０２により対象物理システム１０６を励振させた位置において計測部１０３が計測した入力信号及び出力信号に基づいて自己周波数応答関数を算出する。分析部１０５は、算出された自己周波数応答関数から得られるインパルス応答関数と、分析対象である対象物理システム１０６がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて対象物理システム１０６のシステム同定を行う。これにより、近接固有値を有する系のシステム同定を行うことができる。

Description

システム同定装置、システム同定方法及び記録媒体
　本発明は、システム同定装置、システム同定方法及び記録媒体に関する。
　石油、ガス、水道などのプラントシステムや、産業用ロボット等の物理システムを対象として、ＩｏＴ（Internet of Things）による監視や制御を行う際には、対象とするシステムのモデリング技術が重要である。このモデリング技術として、観測データに基づいて物理システムの数理的モデリングを行うものがある（例えば、特許文献１～４、非特許文献１参照）。
国際公開２０１５／１１８７３７号 国際公開２０１５／０５９９５６号 特開平４－７７７９８号公報 特開平３－２１７９０１号公報
中溝　高好、「信号解析とシステム同定」、コロナ社、ｐ．２２－２４、４９－５３，１２１－１２７、１９８８年
　しかし、モデリングの対象システムが、うなり現象や共鳴が発生するような、近接固有値を有する系である場合、特許文献１～４および非特許文献１に記載の技術では、システム同定が困難な場合があった。なお、近接固有値を有する系は、固有値が縮退した重根を有する系も含むものとする。
　本発明は、近接固有値を有する系のシステム同定を行うことができるシステム同定装置等を提供することを目的としている。
　本発明の第１の態様によれば、システム同定装置は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部と、を備える。
　本発明の第２の態様によれば、システム同定方法は、分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有する。
　本発明の第３の態様によれば、記録媒体は、コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラムを記録する。
　本発明によれば、近接固有値を有する系のシステム同定を行うことができる。
本発明の第１の実施形態によるシステム同定装置の構成を示すブロック図である。 同実施形態によるシステム同定装置の処理を示すフローチャートである。 第２の実施形態によるシステム同定装置の構成を示すブロック図である。 同実施形態によるシステム同定装置の処理を示すフローチャートである。 第３の実施形態によるシステム同定装置の構成を示すブロック図である。 同実施形態による消火栓カプラの状況を示す図である。 同実施形態による仮想の仮想二自由度モデルを示す図である。 同実施形態によるシステム同定結果を示す図である。 システム同定実験の結果を示す図である。 自己回帰モデルを用いた同定方法と本実施形態の比較結果を示す図である。 本発明の実施形態によるシステム同定装置の最小構成を示すブロック図である。
　以下、図面を参照して本発明の一実施の形態を説明する。
　観測データに基づく物理システムの数理的なモデリング手法は、「システム同定問題」と呼ばれている。本問題では大別して、（１）システムの入力および出力の信号が既知の場合と、（２）入力が未知の場合とに分類される。さらには、時間領域信号又は周波数領域信号を用いた技術も知られている。
　時間領域情報を用いた技術では、ＡＲモデル(Autoregressive model：自己回帰モデル)やＭＡモデル(Moving Average model：移動平均モデル)、ＡＲＭＡモデル(Autoregressive moving average model：自己回帰移動平均モデル)、ＡＲＸモデル（Auto-Regressive eXogeneous model）などの多項式モデルが用いられる。多項式モデルでは、周波数領域が平坦化され、異なる固有値が近接している近接固有値を有する系への適用が困難であった。一方で、周波数領域情報を用いた場合には、固有値双方のピーク位置が不明瞭となる場合が多く、曲線適合法を適用することができない。特に、フーリエ変換による周波数領域のサンプル数によっては、ピーク位置を視認することすらできない。
　上記のことから、近接固有値を有する系では固有値が近接するため、周波数領域同定法および時間領域同定法では、システム同定が困難となる。このように、近接固有値を有する系のシステム同定問題は未だ十分に解決されていない問題である。本実施形態のシステム同定装置等は、このような問題を解決し、近接固有値を有する系（固有値が重なる重根を有する系も含む）のシステム同定を行う。
［第１の実施形態］
　図１は、第１の実施形態によるシステム同定装置１の構成を示すブロック図である。システム同定装置１は、設置位置決め部１０１、加振部１０２、計測部１０３、信号収集部１０４及び分析部１０５を有する。対象物理システム１０６は、システム同定装置１による同定対象である。設置位置決め部１０１は、対象物理システム１０６に、加振部１０２及び計測部１０３を設置する。加振部１０２は、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、加振部１０２が設置位置決め部１０１を介して対象物理システム１０６を励振させたときの、対象物理システム１０６への入力信号と出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化する。分析部１０５は、信号収集部１０４により得られたデータを分析し、対象物理システム１０６のシステム同定を行う。
　図２は、システム同定装置１の処理を示すフローチャートである。測定者は、対象物理システム１０６に対し、設置位置決め部１０１を介して加振部１０２及び計測部１０３を設置する（ステップＳ１１０）。このとき、入力位置と出力位置とを一致させるようにする。入力位置とは、加振部１０２が設置される、すなわち、振動の原因が入力される位置である。出力位置とは、計測部１０３が設定される、すなわち、対象物理システム１０６の振動を計測する位置である。
　自己周波数応答関数を計測するため、加振部１０２により、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、対象物理システム１０６への振動の入力信号と対象物理システム１０６からの振動の出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化し、分析部１０５に出力する。分析部１０５は、得られたデータを分析する。具体的には、分析部１０５は、最初に、入力信号と出力信号をそれぞれ高速フーリエ変換（ＦＦＴ：fast fourier transform）する。分析部１０５は、周波数領域で出力信号を入力信号により除算することにより、自己周波数応答関数を求める（ステップＳ１２０）。
　次いで、分析部１０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ１３０）。分析部１０５は、ズーミングした自己周波数応答関数に逆フーリエ変換（Inverse fourier transform）を施すことにより、自己周波数応答関数のインパルス応答関数を求める（ステップＳ１４０）。
　分析部１０５は、次のステップＳ１６０において用いられる初期値及びステップ幅の入力を受ける（ステップＳ１５０）。分析部１０５は、ステップＳ１４０で得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ１６０）。分析部１０５は、多変数ニュートン法の実行時に、ステップＳ１５０において入力された初期値及びステップ幅を用いる。
　分析部１０５は、解が収束しないと判定した場合（ステップＳ１７０：ＮＯ）、ステップＳ１５０に戻って、新たに使用する初期値及びステップ幅の入力を受け、ステップＳ１６０の処理を行う。分析部１０５は、解が収束したと判定した場合（ステップＳ１７０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム１０６の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ１８０）。
　本実施形態によれば、近接固有値を有する系のシステム同定が可能である。
［第２の実施形態］
　本実施形態は、第１の実施形態を管路のシステム同定に適用する。
　図３は、第２の実施形態によるシステム同定装置３の構成を示すブロック図である。システム同定装置３は、設置位置決め部３０１、加振部３０２、計測部３０３、信号収集部３０４、分析部３０５、記憶部３０６及び初期値設定部３０７を有する。対象物理システム３０８は、管路であり、システム同定装置３による同定対象である。
　設置位置決め部３０１、加振部３０２、計測部３０３及び信号収集部３０４はそれぞれ、第１の実施形態の設置位置決め部１０１、加振部１０２、計測部１０３及び信号収集部１０４と同様の機能を有する。記憶部３０６は、管路管理台帳データを記憶する。管路管理台帳データは、対象物理システム３０８の物理データである口径、材種、管厚の情報を含む。初期値設定部３０７は、記憶部３０６から読み出した口径、材種、管厚のデータに基づいて、分析部３０５において用いられる多変数ニュートン法のパラメータの初期値を算出する。分析部３０５は、初期値設定部３０７が算出した初期値を多変数ニュートン法において用いる点以外は、第１の実施形態の分析部１０５と同様の機能を有する。
　図４は、システム同定装置３の処理を示すフローチャートである。測定者は、対象物理システム３０８に対し、設置位置決め部３０１を介して加振部３０２及び計測部３０３を設置する（ステップＳ３１０）。加振部３０２により、設置位置決め部３０１を介して対象物理システム３０８を励振させる。計測部３０３は、加振位置における対象物理システム３０８への入力信号と対象物理システム３０８からの出力信号を検知する。信号収集部３０４は、計測部３０３が検知した入力信号及び出力信号をデータ化する。分析部３０５は、入力信号と出力信号を用いて自己周波数応答関数を求める（ステップＳ３２０）。分析部３０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ３３０）。分析部３０５は、ズーミング部分を逆フーリエ変換し、自己周波数応答関数のインパルス応答関数を求める（ステップＳ３４０）。
　初期値設定部３０７は、記憶部３０６から対象物理システム３０８の口径、材種、管厚のデータを読み出す（ステップＳ３５０）。初期値設定部３０７は、読み出したデータに基づいて多変数ニュートン法のパラメータの初期値を算出する（ステップＳ３６０）。分析部３０５は、ステップ幅の入力を受ける（ステップＳ３７０）。分析部３０５は、ステップＳ３４０において得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ３８０）。多変数ニュートン法においては、ステップＳ３６０において算出された初期値と、ステップＳ３７０において入力されたステップ幅とを用いる。
　分析部３０５は、解が収束しないと判定した場合（ステップＳ３９０：ＮＯ）、ステップＳ３７０に戻って、新たに使用するステップ幅の入力を受け、ステップＳ３８０の処理を行う。分析部３０５は、解が収束したと判定した場合（ステップＳ３９０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム３０８の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ４００）。
　本実施形態によれば、管路のシステム同定が可能である。
［第３の実施形態］
　ここでは、第１の実施形態を水道管路のシステム同定に適用する。
　図５は、本実施形態のシステム同定装置５の構成を示すブロック図である。システム同定装置５は、消火栓カプラ５０１、ハンマー５０２、センサ５０３、データロガー５０４及び同定処理部５０５を有する。管路５０６は、システム同定装置５によるシステム同定対象の水道管路である。ハンマー５０２として、力センサを内蔵したインパルスハンマー、市販ハンマーに加速度ピックアップを設置したもの、電磁式加振器などが例示できる。センサ５０３として、加速度ピックアップ、レーザードップラー型速度計、レーザー変位計、接触式変位計などが例示できる。同定処理部５０５は、例えば、プロセッサ、メモリ及びＨＤＤ（Hard disk drive）により実現される。プロセッサは、処理をコンピュータに実行させるための同定処理プログラムをＨＤＤから読み出して実行することにより同定処理部５０５として動作する。
　図６は、管路５０６に設置された消火栓カプラ５０１の状況を示す図である。計測者は、管路５０６に消火栓カプラ５０１及びセンサ５０３を設置する（図２のステップＳ１１０）。計測者は、ハンマー５０２により図６に示す消火栓カプラ５０１を打撃し、管路５０６を励振する。センサ５０３は、打撃位置（Tapping Point）と同じ測定位置（Measurements point for vibration response）における励振後の出力信号を検出する。データロガー５０４は、ハンマー５０２の入力信号及びセンサ５０３の出力信号を収集する。データロガー５０４は、収集した入力信号及び出力信号をデータ化して同定処理部５０５に出力する。
　同定処理部５０５は、入力信号と出力信号それぞれに高速フーリエ変換（ＦＦＴ）処理を行う。ＦＦＴにより得られた入力信号のスペクトルをＸ（ω）、出力信号のスペクトルをＹ（ω）とそれぞれ表す。ωは周波数である。同定処理部５０５は、各周波数領域のスペクトルを除し、自己周波数応答関数Ｌ（ω）＝Ｙ（ω）／Ｘ（ω）を求める（図２のステップＳ１２０）。ここで、自己周波数応答関数の算出にあたり、Ｌ（ω）＝（Ｙ（ω）・Ｘ＊（ω））／（Ｘ（ω）・Ｘ＊（ω））をＨ１推定、Ｌ（ω）＝（Ｙ（ω）・Ｙ＊（ω））／（Ｘ（ω）・Ｙ＊（ω））をＨ２推定と呼ぶが、どちらを用いても良い。なお、Ｘ＊（ω）はＸ（ω）の複素共役、Ｙ＊（ω）はＹ（ω）の複素共役である。
　同定処理部５０５は、得られた自己周波数応答関数Ｌ（ω）において、着目する近接固有値が存在する周波数帯域のみをズーミングする（図２のステップＳ１３０）。着目する近接固有値が存在する周波数帯域にはピークが現れる。そこで、同定処理部５０５は、自己周波数応答関数Ｌ（ω）におけるピークを検出することにより、着目する近接固有値が存在する周波数帯域を特定する。あるいは、同定処理部５０５は、自己周波数応答関数Ｌ（ω）をシステム同定装置５が備える表示装置に表示し、表示を確認したユーザが、ピークが現れる周波数帯域を入力してもよい。同定処理部５０５は、特定された周波数帯域の自己周波数応答関数Ｌ（ω）を抽出し、閾値よりも小さい値については０に置き換えるズーミングを行う。同定処理部５０５は、ズーミングした結果を逆ＦＦＴすることにより、インパルス応答関数ｇｅ（ｔ）を得る。なお、ｔは時間を表す（図２のステップＳ１４０）。
　ここで、システム同定用モデルとして仮想二自由度モデルを仮想する。なお、仮想二自由度モデルは、対称型二自由度ばね質量系とも呼ばれる。
図７は、仮想二自由度モデルを示す図である。仮想二自由度モデルは、等しい質量、ばね定数、及び、減衰係数からなる一自由度ばね質量系が、ばね及びダッシュポットによって接続された系である。Ｍは質量、Ｋはばね定数、Ｃは減衰係数、Ｆは外力ベクトル、ｘ１、ｘ２は変位ベクトルである。また、ΔＫはばね定数の変化、ΔＣは減衰係数の変化を示す。仮想二自由度モデルは、運動方程式を表す質量行列、剛性行列、減衰行列が対称行列となるシステムであり、固有ベクトルが対称となる性質を有しており、近接固有値を有する系のシステム同定に好適となる。
　図７に示す仮想二自由度モデルの自己周波数応答関数Ｌ１１を求めると、以下の式（１）となる。ｓは複素数を表す。ｓはｓ＝ｉ・ωで表される（ただし、ωは周波数、ｉは虚数単位）。
Figure JPOXMLDOC01-appb-M000001
　式（１）をラブラス逆変換し、仮想二自由度モデルのインパルス応答関数ｇ１１を求めると、式（２）となる。なお、ｓ＝ｉ・ωを式（１）に代入した場合、これを逆フーリエ変換すると、式（２）となる。
Figure JPOXMLDOC01-appb-M000002
　δ（ｔ）はディラックのデルタ関数である。第１項の（１／Ｍ）・δ（ｔ）は、無視できる程度に他の項よりも小さな値である。
　なお、ここで、各パラメータωｄ１、ωｄ２、α、β、γを式（３）とした。
Figure JPOXMLDOC01-appb-M000003
　これより、未知パラメータは計５つとなる。実験値のインパルス応答関数ｇｅ（ｔ）と式（２）の差の二乗和Ｊを目的関数とした多変数ニュートン法によって、パラメータ推定を行う更新式を求める。目的関数Ｊを式（４）に示し、更新式を式（５）に示す。ｉはサンプリングの番号、ｔｉはｉ番目のサンプリングを行った時間を表す。
Figure JPOXMLDOC01-appb-M000004
Figure JPOXMLDOC01-appb-M000005
　ここで、λはステップ調整用パラメータであり、パラメータ推定アルゴリズムが発散する場合に０．００１～０．１の間で調整すると、収束性が改善される。特に、管路系のシステム同定では、λを０．０１程度に設定することが好ましい。ｇ＾１１は、現在のパラメータ値を用いて式（２）を算出したものである。図２のステップＳ１５０においては、各パラメータ（ωｄ１、ωｄ２、α、β、γ）の初期値、及び、ステップ幅となるλの値が入力される。なお、同定処理部５０５は、初期値の算出に用いられる情報の入力を受け、入力された情報に基づいて初期値を算出してもよい。
　同定処理部５０５は、ステップＳ１４０により得られたインパルス応答関数ｇｅ（ｔｉ）と、現在の各パラメータの値を用いて式（２）により算出したｇ＾１１（ｔｉ）とに基づいて、式（４）により二乗和Ｊを算出する。同定処理部５０５は、二乗和Ｊが閾値以下となるように、式（５）により各パラメータの値を更新する（図２のステップＳ１６０）。
　同定処理部５０５は、パラメータの更新を繰り返し、Ｊの値やその値の変化によって収束したか否かを判断する。同定処理部５０５は、収束しないと判断した場合（図２のステップＳ１７０：ＮＯ）、新たな初期値及びλの入力を受ける（ステップＳ１５０）。同定処理部５０５は、収束したと判断した場合、そのときのパラメータ値を用い、式（５）の関係に基づいて、質量Ｍ、ばね定数Ｋ及び減衰係数Ｃを算出する（ステップＳ１８０）。
　数値実験により、本実施形態によるシステム同定法の動作検証をおこなった。口径１００ｍｍのダクタイル鋳鉄管で構成された管路を模擬するため、真値がＭ＝１３．３０７０ｋｇ、Ｋ＝２．８９９４×１０９Ｎ／ｍ、Ｃ＝１０００Ｎｓ／ｍ、ΔＫ＝５．７９８９×１０７Ｎ／ｍ、ΔＣ＝６６６．６６６７Ｎｓ／ｍとした。本条件よりサンプリング周波数５０ｋＨｚで２０ｍｓのインパルス応答関数を生成し、さらに平均０、分散１の正規性白色雑音をインパルス応答関数に加算し、試験用のテストデータとした。初期値は各真値に０．９５を乗じた値を用い、ステップ調整用パラメータは０．０１を用いた。
　図８は、上記の条件によるシステム同定結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。同図には、同定したパラメータを用いて計算した自己周波数応答関数（Identified）と、真値の自己周波数応答関数（Experiment）をそれぞれ示している。同図によれば、真値と、同定結果の両者の良い一致が確認できる。なお、推定値は、Ｍ＝１３．３０７３ｋｇ、Ｋ＝２．８９８０×１０９Ｎ／ｍ、Ｃ＝９９９．９３１２Ｎｓ／ｍ、ΔＫ＝６．０６５２×１０７Ｎ／ｍ、ΔＣ＝６６６．６７３０Ｎｓ／ｍであった。以上により、本実施形態による同定アルゴリズムの動作が確認できる。
　なお、第２の実施形態のように管路のシステム同定を行う場合、パラメータωｄ１、ωｄ２、α、β、γの初期値は以下のように算出される。これらの初期値を決める主たるパラメータである質量Ｍ及びばね定数Ｋは、以下の式（６）を用いて算出される。
Figure JPOXMLDOC01-appb-M000006
　ここで、Ｒは半径、Ａは断面積である。管長をＬ、管厚をｈとした場合、Ａ＝ｈＬの関係がある。半径Ｒ、管厚ｈは管路管理台帳データから読み出した口径、管厚から求められる。管長Ｌは、管路管理台帳データから読み出してもよく、ユーザが入力してもよい。Ｅは管の弾性係数、Ｉは断面二次モーメントであり、Ｉ＝Ｌ・ｈ３／１２である。ρは、管の密度である。弾性係数Ｅ及び管の密度ρは、管路管理台帳データから読み出した材種に応じた値としてもよく、予め決められた値としてもよい。ΔＫはばね定数Ｋの１／１００、ΔＣは減衰係数Ｃの１／３程度が望ましい。初期値の算出に用いる減衰係数Ｃ、減衰係数の変化ΔＣは、管路管理台帳データから読み出した口径、管厚、材種のうち一以上に応じた値としてもよく、予め決められた値としてもよい。初期値設定部３０７は、これらの値を用いて式（３）により各パラメータωｄ１、ωｄ２、α、β、γの初期値を算出する。
以上が初期値設定部の具体的な表式である。
［実験］
　供用後の水道管を試験用管路に設置し、通水環境下でシステム同定実験を実施した。供試管路は口径１００ｍｍ、肉厚１０ｍｍの普通鋳鉄管で供用済みの管を用いた。管路の上部には消火栓を設置し、カプラ上に加速度センサを設置した。
　図９は、システム同定実験の結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。また図中の丸のプロットは自己周波数応答関数の実験値（Experiment）を表し、実線は同定結果（Identified）を表している。同図に示すように、実験値と推定値のピーク位置やスペクトル形状について、良好な一致が確認された。なお、得られた推定値はＭ＝１４．８４４６ｋｇ、Ｋ＝３．４３４６×１０９Ｎ／ｍ、Ｃ＝９０３．５４９１Ｎｓ／ｍ、ΔＫ＝１．２８７７×１０８Ｎ／ｍ、ΔＣ＝９０３．５４９１Ｎｓ／ｍであった。
　関連技術と比較した優位性を示すため、特許文献１で用いられているＡＲモデルを用いた同定を行った。図１０は、ＡＲモデルを用いた同定方法と本実施形態による同定方法との比較結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。
　同図における丸のプロットは自己周波数応答関数の実験値（Experiment）を表し、符号Ｌ１は本実施形態の同定結果（Identified）を表している。また、符号Ｌ２は、ＡＲ法による同定結果（AR）を示す。計算条件はＡＲモデル次数が１００次であり、本実験のインパルス応答関数を同定入力として用いている。これは、本実施形態のシステム同定方法の評価信号と全く同一である。ＡＲ法による同定結果は、実験値との大きな乖離が確認され、同定困難であることが確認できる。これは、インパルス応答関数がうなり波形をともなっており、時間領域同定法で用いる多項式モデルでは、本波形の特徴を記述するのに限界が生じていることを示している。以上の検討により、本実施形態が優れた効果を奏するという結果が示された。
　図１１は、本発明の実施形態によるシステム同定装置の最小構成を示すブロック図である。同図に示す最小構成のシステム同定装置１ａは、上述した第１の実施形態の分析部１０５を少なくとも備えればよい。分析部１０５は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出する。そして、分析部１０５は、算出された自己周波数応答関数から得られるインパルス応答関数と、分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて、分析対象のシステム同定を行う。
　本実施形態によれば、仮想二自由度モデルの自己周波数応答関数に対応するインパルス応答関数を用いて、近接固有値を有する系のシステム同定を行うことができる。
　上述した実施形態におけるシステム同定装置１、１ａ、３、５は、バスで接続されたＣＰＵ（Central Processing Unit）やメモリや補助記憶装置などを備え、システム同定プログラムを実行することによって上述した実施形態におけるシステム同定装置１、１ａ、３、５の一部の機能を実現する。なお、システム同定装置１、１ａ、３、５の一部の機能は、ＡＳＩＣ（Application Specific Integrated Circuit）やＰＬＤ（Programmable Logic Device）やＦＰＧＡ（Field Programmable Gate Array）等のハードウェアを用いて実現されても良い。システム同定プログラムは、コンピュータ読み取り可能な記録媒体に記録されても良い。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ＲＯＭ（Read Only Memory）、ＣＤ－ＲＯＭ（Compact Disc Read only memory）等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。システム同定プログラムは、電気通信回線を介して送信されても良い。
　上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、を備えるシステム同定装置。
（付記２）前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、付記１に記載のシステム同定装置。
（付記３）前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、付記２に記載のシステム同定装置。
（付記４）前記分析対象は、管路である、付記１に記載のシステム同定装置。
（付記５）前記分析対象を励振させる加振部と、前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、付記１に記載のシステム同定装置。
（付記６）前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、付記５に記載のシステム同定装置。
（付記７）前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、付記５に記載のシステム同定装置。
（付記８）前記設置位置決め部は、消火栓カプラである、付記５に記載のシステム同定装置。
（付記９）前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、付記１に記載のシステム同定装置。
（付記１０）分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有するシステム同定方法。
（付記１１）コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラム。
　以上、実施形態（及び実施例）を参照して本願発明を説明したが、本願発明は上記実施形態（及び実施例）に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
　この出願は、２０１８年２月２１日に出願された日本出願特願２０１８－０２９２１８を基礎とする優先権を主張し、その開示の全てをここに取り込む。
　近接固有値を有する系のシステム同定に適用可能である。よってその工業的価値は高い。
１、１ａ、３、５…システム同定装置
１０１、３０１…設置位置決め部
１０２、３０２…加振部
１０３、３０３…計測部
１０４、３０４…信号収集部
１０５、３０５…分析部
１０６、３０８…対象物理システム
３０６…記憶部
３０７…初期値設定部
５０１…消火栓カプラ
５０２…ハンマー
５０３…センサ
５０４…データロガー
５０５…同定処理部
５０６…管路

Claims (11)

  1. 　分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、
    　を備えるシステム同定装置。
  2. 　前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、
    　請求項１に記載のシステム同定装置。
  3. 　前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、
    　請求項２に記載のシステム同定装置。
  4. 　前記分析対象は、管路である、
    　請求項１に記載のシステム同定装置。
  5. 　前記分析対象を励振させる加振部と、
    　前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、
    　前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、
    　請求項１に記載のシステム同定装置。
  6. 　前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、
    　請求項５に記載のシステム同定装置。
  7. 　前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、
    　請求項５に記載のシステム同定装置。
  8. 　前記設置位置決め部は、消火栓カプラである、
    　請求項５に記載のシステム同定装置。
  9. 前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、請求項１に記載のシステム同定装置。
  10. 　分析対象を励振させる加振ステップと、
    　前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、
    　前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、
    　を有するシステム同定方法。
  11. 　コンピュータに、
    　分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、
    　を実行させるプログラムを記録する記録媒体。
PCT/JP2019/005805 2018-02-21 2019-02-18 システム同定装置、システム同定方法及び記録媒体 WO2019163701A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
EP19756575.7A EP3757704A4 (en) 2018-02-21 2019-02-18 SYSTEMS IDENTIFICATION DEVICE, SYSTEM IDENTIFICATION PROCESS AND RECORDING MEDIA
JP2020501746A JP6981526B2 (ja) 2018-02-21 2019-02-18 システム同定装置、システム同定方法及びコンピュータプログラム
US16/970,754 US20210010980A1 (en) 2018-02-21 2019-02-18 System identification device, system identification method, and recording medium

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018029218 2018-02-21
JP2018-029218 2018-02-21

Publications (1)

Publication Number Publication Date
WO2019163701A1 true WO2019163701A1 (ja) 2019-08-29

Family

ID=67688228

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/005805 WO2019163701A1 (ja) 2018-02-21 2019-02-18 システム同定装置、システム同定方法及び記録媒体

Country Status (4)

Country Link
US (1) US20210010980A1 (ja)
EP (1) EP3757704A4 (ja)
JP (1) JP6981526B2 (ja)
WO (1) WO2019163701A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111830062A (zh) * 2020-07-22 2020-10-27 深圳市赛龙自动化科技有限公司 一种陶瓷瓶瑕疵检测设备及其检测方法
CN112765863A (zh) * 2021-02-04 2021-05-07 上海交通大学 基于位姿依赖特性和交叉耦合项的机器人刀尖频响预测方法和系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH03217901A (ja) 1990-01-24 1991-09-25 Toshiba Corp システムの同定装置
JPH0477798A (ja) 1990-07-19 1992-03-11 Oki Electric Ind Co Ltd 周波数包絡線成分の特徴量抽出方法
WO2015059956A1 (ja) 2013-10-23 2015-04-30 日本電気株式会社 構造物診断装置、構造物診断方法、及びプログラム
WO2015118737A1 (ja) 2014-02-07 2015-08-13 三菱電機株式会社 システム同定装置
WO2016114136A1 (ja) * 2015-01-14 2016-07-21 日本電気株式会社 配管検査システム、配管検査装置、配管検査方法、及び、記録媒体
JP2018029218A (ja) 2016-08-15 2018-02-22 ＡＰＲＥＳＩＡ Ｓｙｓｔｅｍｓ株式会社 アンテナ装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050072234A1 (en) * 2003-05-20 2005-04-07 Weidong Zhu System and method for detecting structural damage
EP2682729A1 (en) * 2012-07-05 2014-01-08 Vrije Universiteit Brussel Method for determining modal parameters

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH03217901A (ja) 1990-01-24 1991-09-25 Toshiba Corp システムの同定装置
JPH0477798A (ja) 1990-07-19 1992-03-11 Oki Electric Ind Co Ltd 周波数包絡線成分の特徴量抽出方法
WO2015059956A1 (ja) 2013-10-23 2015-04-30 日本電気株式会社 構造物診断装置、構造物診断方法、及びプログラム
WO2015118737A1 (ja) 2014-02-07 2015-08-13 三菱電機株式会社 システム同定装置
WO2016114136A1 (ja) * 2015-01-14 2016-07-21 日本電気株式会社 配管検査システム、配管検査装置、配管検査方法、及び、記録媒体
JP2018029218A (ja) 2016-08-15 2018-02-22 ＡＰＲＥＳＩＡ Ｓｙｓｔｅｍｓ株式会社 アンテナ装置

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
KOMATSU, MASATAKA ET AL.: "Vibration characteristics estimation for the structure with closely existing eigenvalues by realization theory", JOURNAL OF STRUCTURAL ENGINEERING, JAPAN SOCIETY OF CIVIL ENGINEERS, vol. 59, pages 340 - 352, ISSN: 1881-820X *
NAKAMIZO TAKAYOSHI: "Signal Analysis and System Identification", 1988, CORONA PUBLISHING, pages: 22 - 24,49-53,121-127
NAKAMIZO, TAKAYOSHI ET AL.: "Signal Analysis and System Identification", Tokyo : CORONA PUBUSHING, pages 22 - 24, ISBN: 4-339-03081-3 *
See also references of EP3757704A4

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111830062A (zh) * 2020-07-22 2020-10-27 深圳市赛龙自动化科技有限公司 一种陶瓷瓶瑕疵检测设备及其检测方法
CN112765863A (zh) * 2021-02-04 2021-05-07 上海交通大学 基于位姿依赖特性和交叉耦合项的机器人刀尖频响预测方法和系统
CN112765863B (zh) * 2021-02-04 2022-06-10 上海交通大学 基于位姿依赖特性和交叉耦合项的机器人刀尖频响预测方法和系统

Also Published As

Publication number Publication date
US20210010980A1 (en) 2021-01-14
JPWO2019163701A1 (ja) 2021-02-04
EP3757704A1 (en) 2020-12-30
JP6981526B2 (ja) 2021-12-15
EP3757704A4 (en) 2021-04-07

Similar Documents

Publication Publication Date Title
EP2904368B1 (en) Turbine blade fatigue life analysis using non-contact measurement and dynamical response reconstruction techniques
JP4954729B2 (ja) コンクリート杭の健全性評価支援装置、健全性評価支援方法及び健全性評価支援プログラム
Law et al. Crack identification in beam from dynamic responses
Lam et al. Experimental characterization of multiple cracks in a cantilever beam utilizing transient vibration data following a probabilistic approach
WO2019163701A1 (ja) システム同定装置、システム同定方法及び記録媒体
JP2004069598A (ja) 構造物の損傷推定システムおよびプログラム
JP2010261754A (ja) 常時微動計測に基づく建物の健全性診断法、診断装置及び診断プログラム
Ebrahimzadeh Hassanabadi et al. A Bayesian smoothing for input‐state estimation of structural systems
Debut et al. Identification of the nonlinear excitation force acting on a bowed string using the dynamical responses at remote locations
Smail et al. ARMA models for modal analysis: effect of model orders and sampling frequency
JP6825714B2 (ja) 振動判定装置、振動判定方法及びプログラム
JP2005083975A (ja) 構造性能指標推定装置、及び構造物の構造性能リアルタイムモニタリング方法
JP7004005B2 (ja) 推定装置、推定方法及びコンピュータプログラム
JP2015087172A (ja) 構造物診断装置、構造物診断方法、及びプログラム
Souza et al. Impact of damping models in damage identification
Zhu et al. A two-step approach for structural damage localization and quantification using static and dynamic response data
US20220137003A1 (en) Structure diagnosis apparatus, structure diagnosis method, and computer-readable recording medium
JP2018200217A (ja) 打音検査装置及び打音検査方法
KR101557270B1 (ko) 스마트 구조물의 유지 관리를 위한 est 기반 단일 계측 시스템
Rossing Modal analysis
JP3550296B2 (ja) 構造物の張力および曲げ剛性の測定方法
Manolis et al. Experimental evaluation of damping in beams using the acceleration generalized coordinates: A comparison of the FDD and PCA methods
Rovšček et al. Operational mode-shape normalisation with a structural modification for small and light structures
Sivasubramanian et al. Multiple damage identification in beams using continuous wavelet transforms
US11624687B2 (en) Apparatus and method for detecting microcrack using orthogonality analysis of mode shape vector and principal plane in resonance point

Legal Events

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

Ref document number: 19756575

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2020501746

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2019756575

Country of ref document: EP

Effective date: 20200921