WO2018229898A1 - 状態推定装置 - Google Patents

状態推定装置 Download PDF

Info

Publication number
WO2018229898A1
WO2018229898A1 PCT/JP2017/021964 JP2017021964W WO2018229898A1 WO 2018229898 A1 WO2018229898 A1 WO 2018229898A1 JP 2017021964 W JP2017021964 W JP 2017021964W WO 2018229898 A1 WO2018229898 A1 WO 2018229898A1
Authority
WO
WIPO (PCT)
Prior art keywords
state
equation
time
value
hat
Prior art date
Application number
PCT/JP2017/021964
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 JP2019524625A priority Critical patent/JP6768155B2/ja
Priority to PCT/JP2017/021964 priority patent/WO2018229898A1/ja
Priority to EP17913138.8A priority patent/EP3640822A4/en
Priority to US16/606,770 priority patent/US20200132775A1/en
Publication of WO2018229898A1 publication Critical patent/WO2018229898A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/36—Arrangements for testing, measuring or monitoring the electrical condition of accumulators or electric batteries, e.g. capacity or state of charge [SoC]
    • G01R31/367—Software therefor, e.g. for battery testing using modelling or look-up tables
    • 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
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/048—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators using a predictor
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations

Definitions

  • the present invention relates to a state estimation device that estimates a state or parameter of a target system.
  • Non-Patent Document 1 discloses an apparatus based on an extended Kalman filter (EKF).
  • the extended Kalman filter estimates the internal state of the system based on the output of the target system and a nonlinear discrete-time system that models the target system.
  • the state and the parameter are estimated by converting to an expanded system including the unknown parameter as a state and applying an extended Kalman filter.
  • Such an extended Kalman filter applies the Kalman filter to a nonlinear system by linearizing the system by first-order Taylor approximation, and requires a Jacobian of a state equation or an observation equation.
  • the target system is a non-linear discrete time system, it is necessary to discretize in advance if the state equation is given in continuous time.
  • Patent Document 1 discloses a technique for estimating a state and an unknown parameter by applying an extended Kalman filter by discretizing in advance by the Euler method when a target system is given by a nonlinear continuous time system.
  • the extended Kalman filter when the target system is a nonlinear discrete-time system and the state equation is given in continuous time, it needs to be discretized. In the state estimation by the extended Kalman filter, it is necessary to obtain the Jacobian of the state equation in the process.
  • the present invention has been made in view of the above, and an object of the present invention is to obtain a state estimation device capable of estimating a state with high accuracy by an extended Kalman filter even if the state equation of the target system is a nonlinear continuous-time function. To do.
  • the present invention provides a state in which the state of the target system is estimated by an extended Kalman filter using the output of the target system and a nonlinear model obtained by modeling the target system as inputs.
  • the nonlinear model includes a nonlinear continuous-time state equation, obtains a state estimated value obtained by estimating the state of the target system at a certain time based on the observation equation of the nonlinear model, and an error covariance of the state estimated value State and error estimation unit for obtaining an estimated value of a matrix, and state equation discrete for obtaining a discrete-time state equation by discretizing the nonlinear continuous-time state equation by a quadratic or higher integration method based on the state estimated value at the time And an approximate value of the Jacobian of the time by a difference method using the discrete time state equation of the time and the discrete time equation.
  • a state equation linearization unit a state that predicts the error covariance matrix after a minute time of the time from the approximate value of the Jacobian, and a state that predicts the state after the minute time of the time from the discrete time state equation And an error prediction unit.
  • the state estimation device capable of estimating the state with high accuracy by the extended Kalman filter.
  • FIG. Block diagram showing an outline of the overall configuration of the state estimation apparatus according to Embodiment 1 Block diagram showing details of the state and error prediction unit shown in FIG.
  • FIG. Block diagram showing an outline of the overall configuration of the state estimation apparatus according to the second embodiment The flowchart which shows operation
  • FIG. Block diagram showing an outline of the overall configuration of the state estimation apparatus according to the third embodiment Block diagram showing an outline of the overall configuration of a state estimation apparatus according to Embodiment 4
  • each parameter is a vector quantity except for the time t and the scalar value ⁇ , but is not particularly distinguished in the specification, and the vector quantity is shown in bold in the mathematical expression.
  • the time differentiation is expressed as (d / dt) x in the specification, but a dot is added on the character in the formula.
  • the predicted value and the estimated value are represented by (hat) in the specification, but a hat symbol is added on the character in the mathematical expression.
  • the state quantity or equation of the expanded system model is expressed by (tilde) in the specification, but a tilde is added on the character in the mathematical expression.
  • FIG. 1 is a block diagram showing an outline of the overall configuration of the state estimation apparatus according to the present embodiment.
  • the state estimation apparatus 10 illustrated in FIG. 1 includes a state and error estimation unit 1, a state equation discretization unit 2, a state equation linearization unit 3, and a state and error prediction unit 4.
  • the state and error estimator 1 includes the coefficient C of the linear observation equation, the system output y (t k ) at the current time t k , the observation noise covariance matrix R, and the current state prediction calculated at the previous time. Based on the value x k
  • the state equation discretization unit 2 is a discrete discretized by the fourth-order Runge-Kutta method based on the nonlinear continuous-time state equation f (t, x) and the estimated value x k
  • k (hat)) is output.
  • the state equation linearization unit 3 calculates the Jacobian approximation value ⁇ f d (t, x) / ⁇ x T
  • t t from the discrete-time state equation f d (t k , x k
  • k (hat)). k, x x k
  • the subscript T represents a transposed matrix.
  • FIG. 2 is a block diagram showing details of the state and error prediction unit 4 shown in FIG.
  • the state and error prediction unit 4 illustrated in FIG. 2 includes a state prediction unit 5 and an error prediction unit 6.
  • the state prediction unit 5 outputs the discrete-time state equation f d (t k , x k
  • the error prediction unit 6 calculates ⁇ f d (t, x) / ⁇ x T
  • FIG. 3 is a block diagram showing a hardware configuration of the state estimation device 10 shown in FIG.
  • the state estimation apparatus 10 shown in FIG. 3 includes a processor 21 that performs state estimation, a memory 22 that stores programs and data, and a system interface 23 that reads the output of the target system, and inputs initial values and design parameters.
  • the input device 24 is connected to a display device 25 which is an example of an estimation result output unit.
  • the processor 21 executes a program stored in the memory 22 based on the initial value and design parameters set via the input device 24 and the output value of the target system obtained via the system interface 23.
  • the estimated state value of the target system is calculated.
  • the estimated state value output from the processor 21 is displayed on the display device 25.
  • the processor 21 is, for example, a CPU (Central Processing Unit)
  • the memory 22 is, for example, a RAM (Random Access Memory), a ROM (Read Only Memory), or a flash memory.
  • FIG. 4 is a flowchart showing the operation of the state estimation apparatus according to the present embodiment.
  • the nonlinear model of the target system is assumed to be represented by the following equation (1) which is a nonlinear state equation and a linear observation equation.
  • the time t a value obtained by estimating the state of k x k
  • the state and error estimation unit 1 first performs initialization (S1). Specifically, the state and error estimation unit 1 uses, as initial values, the predicted state value x 1
  • the state and the error estimation unit 1 calculates the Kalman gain K k (S2). Specifically, status and error estimation unit 1 at time t k, the predicted value sigma k of the previous time t k-1 to the calculated current error covariance matrix
  • the Kalman gain K k is calculated by the following equation (2).
  • k ⁇ 1 (hat) is an n ⁇ n matrix
  • R is a q ⁇ q matrix
  • K k is an n ⁇ q matrix.
  • the subscript -1 represents an inverse matrix.
  • the state and error estimation unit 1 calculates an error covariance matrix estimation value ⁇ k
  • k (hat) is an n ⁇ n matrix.
  • the state and error estimation unit 1 calculates a state estimated value x k
  • k (hat) is an n ⁇ 1 matrix.
  • the state x (t k + 1 ) at the next time is expressed by the following equation (5) according to the current state x (t k ) and the sum of integrals from time t k to t k + 1 of the state equation: Become.
  • a calculation method using a fourth-order Runge-Kutta method is shown below.
  • ⁇ f (t k , x (t k )) is calculated from the current state x (t k ).
  • ⁇ f (t k + ( ⁇ T / 2), ⁇ ) is calculated using the calculation result of ⁇ .
  • ⁇ x (t k ) + ( ⁇ T / 2) ⁇ is calculated using the calculation result of ⁇ .
  • ⁇ f (t k + ( ⁇ T / 2), ⁇ ) is calculated using the calculation result of ⁇ .
  • ⁇ x (t k ) + ⁇ T ⁇ is calculated using the calculation result of ⁇ .
  • ⁇ f (t k + 1 , ⁇ ) is calculated using the calculation result of ⁇ .
  • the state equation discretization unit 2 substitutes x k
  • the state prediction unit 5 calculates a state prediction value x k + 1
  • the state equation linearization unit 3 obtains the Jacobian approximate value F k by a difference method based on the discrete time state equation f d and a predetermined small scalar value ⁇ .
  • F k is an n ⁇ n matrix.
  • the state equation linearization unit 3 calculates x i (hat) (S8). Specifically, the state equation linearization unit 3 expresses a value obtained by multiplying only the i-th component of the state by (1 + ⁇ ) from the current state estimated value x k
  • the state equation linearization unit 3 performs calculation by a second-order or higher numerical integration method (S9), and calculates the i-th column of the Jacobian approximate value by the difference method (S10). Specifically, the state equation linearization unit 3 obtains the value of the discrete-time state equation based on x i (hat) by the above-described numerical integration method, and in the current state estimated value x k
  • the error prediction unit 6 calculates an error covariance matrix predicted value ⁇ k + 1
  • Q is an n ⁇ n matrix.
  • the state estimating device 10 advances the time from t k to t k + 1, by repeating to S11 Step S2 described above, and outputs the state estimate of each time.
  • the calculation processing in the above steps is executed by the processor 21 and the Kalman gain K k , estimated value x k
  • the state estimation device even when the state equation of the target system is given by a nonlinear continuous time function, the approximate value of the Jacobian is calculated based on the state estimated value at each time.
  • the discretization error in the process of state estimation by the extended Kalman filter can be reduced, and the accuracy of state estimation can be improved.
  • a state estimation device using an unscented Kalman filter in addition to an extended Kalman filter is known. Since the state estimation apparatus according to the present invention is based on the extended Kalman filter, it requires fewer design parameters than the unscented Kalman filter and can reduce the amount of calculation.
  • the fourth-order Runge-Kutta method is used for approximation of the above equation (5).
  • the present invention is not limited to this, and the two methods such as the Hein method are used. What is necessary is to use an integration method with the accuracy of the next or higher accuracy.
  • the input device It may be possible to adjust from the outside via 24.
  • the fine scalar value ⁇ is adjusted from the outside, the error due to the Jacobian approximation can be reduced by setting the minimum value within a range in which an abnormality such as divergence of the estimation result does not occur.
  • ⁇ may be a matrix of vector n ⁇ 1, and may be obtained by using a minute value ⁇ (i) that is partially or entirely different for each state as in the following equation (10). Good.
  • FIG. FIG. 5 is a block diagram showing an outline of the overall configuration of the state estimation apparatus according to the present embodiment.
  • a state estimation device 10A shown in FIG. 5 includes an observation equation linearization unit 7 and a state estimation device 10 shown in FIG. Note that the description of the same configuration and operation as in the first embodiment will be omitted.
  • the system output y (t k ) at each time are input, and the system noise covariance matrix Q and the observation noise covariance matrix R are used as design parameters.
  • the observation equation linearization unit 7 analytically calculates the Jacobian H from the observation equation h (t, x).
  • the state and error estimator 1 includes an observation equation Jacobian H k , system output y (t k ), current state predicted values x k
  • the state equation discretization unit 2, the state equation linearization unit 3, and the state and error prediction unit 4 operate in the same manner as in the first embodiment.
  • FIG. 6 is a flowchart showing the operation of the state estimation apparatus according to the present embodiment.
  • FIG. 4 is used about S1, S2 to S11.
  • the operation of the observation equation linearization unit 7 is added between step S1 and step S2 in the flowchart shown in FIG.
  • the observation equation linearization unit 7 analytically calculates the Jacobian H of the observation equation as shown in the following equation (13) (S21).
  • H is a q ⁇ n matrix.
  • the observation equation linearization unit 7 substitutes the state predicted value x k
  • k ⁇ 1 (hat) is an n ⁇ 1 matrix
  • H k is a q ⁇ n matrix.
  • the observation equation can be approximated as a value obtained by multiplying the state by the Jacobian around the current state predicted value x k
  • the coefficient C of the observation equation in the first embodiment is replaced with the Jacobian H k , and the output nonlinear prediction function y k
  • the calculation processing in the above steps is executed by the processor 21 and the Kalman gain K k , estimated value x k
  • the state estimation device even when the state equation of the target system is a nonlinear continuous-time function and the observation equation is given by a nonlinear function, there is no need for discretization. Analyzing the Jacobian of the observation equation analytically reduces the amount of calculation, and the Jacobian of the state equation obtains an approximate value based on the estimated state value at each time, so that the discriminating of the state equation can be performed at a second or higher order.
  • the integration method can be applied. Thereby, the discretization error in the process of state estimation by the extended Kalman filter can be reduced, and the accuracy of state estimation can be improved.
  • FIG. 7 is a block diagram showing an outline of the overall configuration of the state estimation apparatus according to the present embodiment.
  • the state estimation apparatus 10B shown in FIG. 7 includes an observation equation linearization unit 7, a state and error estimation unit 1, a state equation discretization unit 2, a state equation linearization unit 3, a state prediction unit 5, and an error prediction. Part 6. Note that the description of the same configuration and operation as in the first and second embodiments is omitted.
  • the hardware configuration uses FIG. 3, and the flowchart uses FIG.
  • the observation equation linearization unit 7 analytically calculates the Jacobian H from the observation equation h (t, x, u).
  • the state and error estimation unit 1 includes an observation equation Jacobian H k , system output y (t k ), current state predicted value x k
  • k (hat) of the state and the error covariance matrix are estimated from k ⁇ 1 (hat).
  • the state equation linearization unit 3 uses the discretized time state equation f d (t k , x k
  • k (hat), u u (t k ) are calculated by the difference method.
  • the state prediction unit 5 outputs the discretized time state equation f d (t k , x k
  • the error prediction unit 6 calculates the Jacobian ⁇ f d (t, x, u) / ⁇ x T
  • k (hat), u u (t k ) of the state equation, From the estimated value ⁇ k
  • the initialization of S1 is the same as in the first embodiment.
  • the observation equation linearization unit 7 analytically calculates the Jacobian H of the observation equation as shown in the following equation (17) (S21).
  • H is a q ⁇ n matrix.
  • the observation equation linearization unit 7 substitutes the predicted value x k
  • k ⁇ 1 (hat) is an n ⁇ 1 matrix
  • H k is a q ⁇ n matrix.
  • the state and error estimator 1 calculates the following equation (19) using the predicted value ⁇ k
  • the state and error estimation unit 1 sets the current error covariance matrix estimated value ⁇ k
  • k ⁇ 1 (hat) by the Kalman gain K k and the coefficient H k of the observation equation is calculated as the following equation (20) (S3).
  • k (hat) is an n ⁇ n matrix.
  • the state and error estimation unit 1 determines the current state estimated value x k
  • k ⁇ 1 (hat) predicted from the observation equation based on (t k ) by the Kalman gain is the current state predicted value x k
  • the added value is calculated as in the following equation (21) (S4).
  • k (hat) is an n ⁇ 1 matrix.
  • a calculation method based on a fourth-order Runge-Kutta method considering system input is shown below.
  • ⁇ f (t k , x (t k ), u (t k )) is calculated from the current state x (t k ) and the system input u (t k ).
  • ⁇ f (t k + ( ⁇ T / 2), ⁇ , u (t k )) is calculated using the calculation result of ⁇ .
  • ⁇ x (t k ) + ( ⁇ T / 2) ⁇ is calculated using the calculation result of ⁇ .
  • ⁇ f (t k + ( ⁇ T / 2), ⁇ , u (t k )) is calculated using the calculation result of ⁇ .
  • ⁇ x (t k ) + ⁇ T ⁇ is calculated using the calculation result of ⁇ .
  • ⁇ f (t k + 1 , ⁇ , u (t k )) is calculated using the calculation result of ⁇ .
  • the state equation discretization unit 2 substitutes the current state estimated value x k
  • k (hat), u (t k )) are calculated (S6).
  • the state prediction unit 5 converts the discrete time state equation f d (t k , x k
  • k (hat) is output (S7).
  • the state equation linearization unit 3 generates a Jacobian approximate value F based on a discrete time state equation f d (t k , x k
  • F k is an n ⁇ n matrix.
  • the state equation linearization unit 3 obtains x i (hat) from the current state estimated value x k
  • the state equation linearization unit 3 obtains the value of the discrete time state equation by the above-described numerical integration method based on x i (hat) (S9), and the current state estimated value x k
  • k (i ) (Hat), the Jacobian i-th column (i 1, 2,..., N) is approximated by the following equation (24) (S10).
  • the error prediction unit 6 calculates the following equation (6) based on the Jacobian approximate value F k , the current error covariance matrix estimate ⁇ k
  • Q is an n ⁇ n matrix.
  • the calculation processing in the above steps is executed by the processor 21 and the Kalman gain K k , estimated value x k
  • the state estimation device even when the target system has an input and the state equation of the target system is given by a nonlinear continuous-time function, the state estimation value at each time is obtained.
  • the state estimation value at each time is obtained.
  • the discretization error in the process of state estimation by the extended Kalman filter can be reduced, and the accuracy of state estimation can be improved.
  • the fourth-order Runge-Kutta method is used for approximation of the above equation (22), but the present invention is not limited to this, and the Hoin method or the like is used. An integration method with second or higher accuracy may be used.
  • the observation equation of the target system is nonlinear.
  • the present invention can also be applied when the observation equation is a linear function.
  • the observation equation linearization unit 7 may be omitted, and the state coefficient matrix in the observation equation may be H k .
  • ⁇ may be a matrix of vector n ⁇ 1, and may be obtained using a minute value ⁇ (i) that is partially or entirely different for each state as in the following equation (26). Good.
  • FIG. 8 is a block diagram showing an outline of the overall configuration of the state estimation apparatus according to the present embodiment.
  • a state estimation device 10C illustrated in FIG. 8 includes a model conversion unit 8 and a state estimation device 10B illustrated in FIG. Note that the description of the same configuration and operation as those in Embodiments 1 to 3 is omitted.
  • the hardware configuration uses FIG. 3, and the flowchart uses FIG.
  • the system output y (t k ) at each time is input, the system noise covariance matrix Q, the observation noise covariance matrix R, and the small scalar value ⁇ in the difference method are used as design parameters, and the unknown parameter ⁇ and the state x are set.
  • the model conversion unit 8 converts the continuous time model into a nonlinear continuous time model using the original state x and the unknown parameter ⁇ as a new state vector.
  • the state estimation device 10B receives the converted nonlinear continuous time model as an input and calculates a new state vector estimate.
  • the observation equation linearization unit 7, the state and error estimation unit 1, the state equation discretization unit 2, the state equation linearization unit 3, the state prediction unit 5 and the error prediction unit 6 in the state estimation device 10B are the same as those in the first to third embodiments. Is the same.
  • x (tilde) is an (n + p) ⁇ 1 matrix.
  • the unknown parameter ⁇ it is modeled to be the n + 1: n + p component of f (tilde).
  • the n + 1: n + p component of f (tilde) is set to 0.
  • the subsequent operation is the same as in FIG. 4, and an estimated value of a new vector x (tilde) is calculated by the following equation (30).
  • the calculation process in the above steps is executed by the processor 21 and the Kalman gain K k , the estimated value x k
  • k ⁇ 1 (hat), discrete time state equation f d (tilde), and Jacobian approximate value F k (tilde) of the state equation are stored in the memory 22. In each step, it is read from the memory 22 and used.
  • the state estimation device even when the state equation of the target system is a continuous-time function and includes an unknown parameter, the original state and the unknown parameter are set to a new state. Conversion to a model makes it possible to estimate the state, reduce discretization errors, and perform state estimation and parameter estimation with high accuracy.
  • the continuous-time model including the unknown parameter as an input of the present embodiment may be nonlinear or linear. Even in the linear case, the output of the model conversion unit is a nonlinear model.
  • the unknown parameter may be included in one or both of the state equation and the observation equation, and may be a continuous time model in which the output of the model conversion unit includes a nonlinear continuous time state equation.
  • the continuous time model may or may not include system inputs.
  • the configuration described in the above embodiment shows an example of the contents of the present invention, and can be combined with another known technique, and can be combined with other configurations without departing from the gist of the present invention. It is also possible to omit or change the part.
  • 1 state and error estimation unit 2 state equation discretization unit, 3 state equation linearization unit, 4 state and error prediction unit, 5 state prediction unit, 6 error prediction unit, 7 observation equation linearization unit, 8 model conversion unit, 10, 10A, 10B, 10C State estimation device, 21 processor, 22 memory, 23 system interface, 24 input device, 25 display device.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Automation & Control Theory (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Health & Medical Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

対象システムの出力と、対象システムをモデル化した非線形モデルとを入力として、拡張カルマンフィルタによって対象システムの状態を推定する状態推定装置において、非線形モデルは非線形連続時間状態方程式を含み、非線形モデルの観測方程式に基づき、ある時刻の対象システムの状態を推定した状態推定値を求め、状態推定値の誤差共分散行列の推定値を求める状態及び誤差推定部と、この時刻の状態推定値に基づいて、非線形連続時間状態方程式を二次以上の積分手法で離散化して離散時間状態方程式を求める状態方程式離散化部と、この時刻の離散時間状態方程式を用いて、差分法によってこの時刻のヤコビアンの近似値を求める状態方程式線形化部と、ヤコビアンの近似値からこの時刻の微小時間後の誤差共分散行列を予測し、離散時間状態方程式からこの時刻の微小時間後の状態を予測する状態及び誤差予測部とを備える。

Description

状態推定装置
　本発明は、対象システムの状態又はパラメータを推定する状態推定装置に関する。
　従来の非線形システムの状態推定装置として、例えば、非特許文献１には拡張カルマンフィルタ（ＥＫＦ：Extended　Kalman　Filter）に基づくものが開示されている。拡張カルマンフィルタは、対象システムの出力と、対象システムをモデル化した非線形離散時間システムとに基づいて、システムの内部状態を推定する。また、システムに未知パラメータを含む場合には、未知パラメータを状態として含む拡大系に変換し、拡張カルマンフィルタを適用することで、状態及びパラメータを推定する。このような拡張カルマンフィルタは、一次テイラー近似によりシステムを線形化することで、非線形システムにカルマンフィルタを適用するものであり、状態方程式又は観測方程式のヤコビアンを要する。また、対象システムを非線形離散時間システムとしているため、状態方程式が連続時間で与えられる場合には、予め離散化することを要する。
　例えば、特許文献１には、対象システムが非線形連続時間システムで与えられる場合に、オイラー法により予め離散化することで拡張カルマンフィルタを適用し、状態及び未知パラメータを推定する技術が開示されている。
　このように、拡張カルマンフィルタでは、対象システムを非線形離散時間システムとし、状態方程式が連続時間で与えられる場合には、離散化することを要する。拡張カルマンフィルタによる状態推定では、その過程で、状態方程式のヤコビアンを求めることを要する。
特開２００５－２４８９４６号公報
谷萩隆嗣、「カルマンフィルタと適応信号処理」、コロナ社、２００５年１２月、ｐ．４８－５１
　しかしながら、上記の従来技術では、連続時間状態方程式の離散化に二次以上の離散化手法を用いた場合に、ヤコビアンを求めることが困難となる。従って、解析的にヤコビアンの導出が可能なオイラー法により離散化することを要するが、オイラー法は一次の近似精度であるため、対象が連続時間状態方程式で与えられる場合には、離散化誤差が大きく、精度の低い推定しかできない、という問題があった。
　本発明は、上記に鑑みてなされたものであって、対象システムの状態方程式が非線形連続時間関数であっても、拡張カルマンフィルタにより高精度に状態を推定可能な状態推定装置を得ることを目的とする。
　上述した課題を解決し、目的を達成するために、本発明は、対象システムの出力と、前記対象システムをモデル化した非線形モデルとを入力として、拡張カルマンフィルタによって前記対象システムの状態を推定する状態推定装置において、前記非線形モデルは非線形連続時間状態方程式を含み、前記非線形モデルの観測方程式に基づき、ある時刻の前記対象システムの状態を推定した状態推定値を求め、前記状態推定値の誤差共分散行列の推定値を求める状態及び誤差推定部と、前記時刻の前記状態推定値に基づいて、前記非線形連続時間状態方程式を二次以上の積分手法で離散化して離散時間状態方程式を求める状態方程式離散化部と、前記時刻の前記離散時間状態方程式を用いて、差分法によって前記時刻のヤコビアンの近似値を求める状態方程式線形化部と、前記ヤコビアンの前記近似値から前記時刻の微小時間後の前記誤差共分散行列を予測し、前記離散時間状態方程式から前記時刻の微小時間後の前記状態を予測する状態及び誤差予測部とを備えることを特徴とする。
　本発明によれば、対象システムの状態方程式が非線形連続時間関数であっても、拡張カルマンフィルタにより高精度に状態を推定可能な状態推定装置を得ることができるという効果を奏する。
実施の形態１に係る状態推定装置の全体構成の概略を示すブロック線図 図１に示す状態及び誤差予測部の詳細を示すブロック線図 図１に示す状態推定装置のハードウェア構成を示すブロック図 実施の形態１に係る状態推定装置の動作を示すフローチャート 実施の形態２に係る状態推定装置の全体構成の概略を示すブロック線図 実施の形態２に係る状態推定装置の動作を示すフローチャート 実施の形態３に係る状態推定装置の全体構成の概略を示すブロック線図 実施の形態４に係る状態推定装置の全体構成の概略を示すブロック線図
　以下に、本発明の実施の形態に係る状態推定装置を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。
　なお、以下の説明において、各パラメータは、時刻ｔ及びスカラー値εを除いてベクトル量であるが、明細書中においては特に区別せずに表記し、数式中ではベクトル量を太字で表記する。また、時間微分は、明細書中においては（ｄ／ｄｔ）ｘで表記するが、数式中では当該文字の上にドットを付する。また、予測値及び推定値は、明細書中においては（ハット）で表記するが、数式中では当該文字の上にハット記号を付する。また、拡大系モデルの状態量又は方程式は、明細書中においては（チルダ）で表記するが、数式中では当該文字の上にチルダを付する。
実施の形態１．
　図１は、本実施の形態に係る状態推定装置の全体構成の概略を示すブロック線図である。図１に示す状態推定装置１０は、状態及び誤差推定部１と、状態方程式離散化部２と、状態方程式線形化部３と、状態及び誤差予測部４とを備える。本実施の形態に係る状態推定装置では、対象システムの非線形連続時間状態方程式（ｄ／ｄｔ）ｘ＝ｆ（ｔ，ｘ）と、線形の観測方程式ｙ＝Ｃｘと、各時刻のシステム出力ｙ（ｔｋ）とを入力とし、システム雑音共分散行列Ｑと、観測雑音共分散行列Ｒとを設計パラメータとする。
　状態及び誤差推定部１は、線形の観測方程式の係数Ｃと、現在の時刻ｔｋにおけるシステム出力ｙ（ｔｋ）と、観測雑音共分散行列Ｒと、前の時刻に算出した現在の状態予測値ｘｋ｜ｋ－１（ハット）と、誤差共分散行列予測値Σｋ｜ｋ－１（ハット）とに基づき、現在の状態と誤差共分散行列の推定値ｘｋ｜ｋ（ハット），Σｋ｜ｋ（ハット）とを推定する。
　状態方程式離散化部２は、非線形連続時間状態方程式ｆ（ｔ，ｘ）と、現在の状態の推定値ｘｋ｜ｋ（ハット）とに基づいて、四次のルンゲクッタ法により離散化された離散時間状態方程式ｆｄ（ｔｋ，ｘｋ｜ｋ（ハット））を出力する。
　状態方程式線形化部３は、離散時間状態方程式ｆｄ（ｔｋ，ｘｋ｜ｋ（ハット））から状態方程式のヤコビアンの近似値∂ｆｄ（ｔ，ｘ）／∂ｘＴ｜ｔ＝ｔｋ、ｘ＝ｘｋ｜ｋ（ハット）を差分法により算出する。なお、添字Ｔは転置行列を表す。
　図２は、図１に示す状態及び誤差予測部４の詳細を示すブロック線図である。図２に示す状態及び誤差予測部４は、状態予測部５と誤差予測部６とを備える。
　状態予測部５は、離散時間状態方程式ｆｄ（ｔｋ，ｘｋ｜ｋ（ハット））を次の時刻の状態の予測値ｘｋ＋１｜ｋ（ハット）として出力する。誤差予測部６は、状態方程式のヤコビアンＦｋである∂ｆｄ（ｔ，ｘ）／∂ｘＴ｜ｔ＝ｔｋ、ｘ＝ｘｋ｜ｋ（ハット）と、現在の誤差共分散行列の推定値Σｋ｜ｋ（ハット）と、システム雑音共分散行列Ｑとから、次の時刻の誤差共分散行列の予測値Σｋ＋１｜ｋ（ハット）を算出する。
　図３は、図１に示す状態推定装置１０のハードウェア構成を示すブロック図である。図３に示す状態推定装置１０は、状態推定を実行するプロセッサ２１と、プログラム及びデータを格納するメモリ２２と、対象システムの出力を読み込むシステムインターフェース２３とを備え、初期値及び設計パラメータを入力する入力装置２４と、推定結果の出力手段の一例である表示装置２５とに接続される。プロセッサ２１は、入力装置２４を介して設定された初期値及び設計パラメータと、システムインターフェース２３を介して取得した対象システムの出力値とに基づいて、メモリ２２に格納されたプログラムを実行することで対象システムの状態推定値を算出する。プロセッサ２１から出力された状態推定値は、表示装置２５に表示される。プロセッサ２１は、例えばＣＰＵ（Central　Processing　Unit）であり、メモリ２２は、例えばＲＡＭ（Random　Access　Memory）、ＲＯＭ（Read　Only　Memory）又はフラッシュメモリである。
　図４は、本実施の形態に係る状態推定装置の動作を示すフローチャートである。対象システムの非線形モデルは、非線形の状態方程式及び線形の観測方程式である下記の式（１）で表されるものとする。
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
　上記ステップにおける計算処理は、プロセッサ２１で実行され、各ステップで求めたカルマンゲインＫｋ、推定値ｘｋ｜ｋ（チルダ）（ハット），Σｋ｜ｋ（ハット）、予測値ｘｋ｜ｋ－１（チルダ）（ハット），Σｋ｜ｋ－１（ハット）、離散時間状態方程式ｆｄ（チルダ）及び状態方程式のヤコビアンの近似値Ｆｋ（チルダ）等の算出値はメモリ２２に格納され、各ステップにおいてメモリ２２から読み出されて利用される。
　以上説明したように、本実施の形態に係る状態推定装置によれば、対象システムの状態方程式が連続時間関数で、未知パラメータを含む場合にも、元の状態と未知パラメータとを新しい状態とするモデルに変換することで状態の推定が可能であり、離散化誤差を軽減し、高精度に状態推定及びパラメータ推定が可能である。
　なお、本実施の形態の入力とする未知パラメータを含む連続時間モデルは、非線形でも線形でもよく、線形の場合でも、モデル変換部の出力は非線形モデルとなる。また、未知パラメータは、状態方程式及び観測方程式のいずれか一方又は双方に含まれていてもよく、モデル変換部の出力が非線形連続時間状態方程式を含むような連続時間モデルであればよい。また、連続時間モデルは、システム入力を含んでもよいし、含まなくてもよい。
　以上の実施の形態に示した構成は、本発明の内容の一例を示すものであり、別の公知の技術と組み合わせることも可能であるし、本発明の要旨を逸脱しない範囲で、構成の一部を省略、変更することも可能である。
　１　状態及び誤差推定部、２　状態方程式離散化部、３　状態方程式線形化部、４　状態及び誤差予測部、５　状態予測部、６　誤差予測部、７　観測方程式線形化部、８　モデル変換部、１０，１０Ａ，１０Ｂ，１０Ｃ　状態推定装置、２１　プロセッサ、２２　メモリ、２３　システムインターフェース、２４　入力装置、２５　表示装置。

Claims (6)

  1. 　対象システムの出力と、前記対象システムをモデル化した非線形モデルとを入力として、拡張カルマンフィルタによって前記対象システムの状態を推定する状態推定装置において、
    　前記非線形モデルは非線形連続時間状態方程式を含み、
    　前記非線形モデルの観測方程式に基づき、ある時刻の前記対象システムの状態を推定した状態推定値を求め、前記状態推定値の誤差共分散行列の推定値を求める状態及び誤差推定部と、
    　前記時刻の前記状態推定値に基づいて、前記非線形連続時間状態方程式を二次以上の積分手法で離散化して離散時間状態方程式を求める状態方程式離散化部と、
    　前記時刻の前記離散時間状態方程式を用いて、差分法によって前記時刻のヤコビアンの近似値を求める状態方程式線形化部と、
    　前記ヤコビアンの前記近似値から前記時刻の微小時間後の前記誤差共分散行列を予測し、前記離散時間状態方程式から前記時刻の微小時間後の前記状態を予測する状態及び誤差予測部とを備える状態推定装置。
  2. 　差分法に用いる微小スカラー値を入力とし、
    　前記状態方程式線形化部が、前記離散時間状態方程式から前記微小スカラー値に基づいて、差分法により前記ヤコビアンの前記近似値を求めることを特徴とする請求項１に記載の状態推定装置。
  3. 　差分法に用いる微小値ベクトルを入力とし、
    　前記状態方程式線形化部が、前記離散時間状態方程式から、各状態に応じた状態方程式の偏微分の近似値を、対応する前記微小値ベクトルの各成分に基づいて、差分法により求めることを特徴とする請求項１に記載の状態推定装置。
  4. 　前記非線形モデルが非線形の観測方程式を含む場合には、前記非線形の観測方程式のヤコビアンを解析的に求める観測方程式線形化部を前記状態及び誤差推定部の前段に備えることを特徴とする請求項１から３のいずれか一項に記載の状態推定装置。
  5. 　前記対象システムの各時刻の入力と前記状態推定値とに基づいて前記非線形連続時間状態方程式を二次以上の積分手法で離散化するモデル離散化部を備えることを特徴とする請求項１から４のいずれか一項に記載の状態推定装置。
  6. 　前記対象システムをモデル化した連続時間モデルが未知パラメータを含み、
    　前記未知パラメータと元の状態とを合わせたベクトルを新しい状態とする拡大系非線形連続時間モデルに変換するモデル変換部を備え、
    　入力となる前記非線形モデルを前記モデル変換部で変換された前記拡大系非線形連続時間モデルとすることを特徴とする請求項１から５のいずれか一項に記載の状態推定装置。
PCT/JP2017/021964 2017-06-14 2017-06-14 状態推定装置 WO2018229898A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2019524625A JP6768155B2 (ja) 2017-06-14 2017-06-14 状態推定装置
PCT/JP2017/021964 WO2018229898A1 (ja) 2017-06-14 2017-06-14 状態推定装置
EP17913138.8A EP3640822A4 (en) 2017-06-14 2017-06-14 DEVICE FOR ESTIMATING THE CONDITION
US16/606,770 US20200132775A1 (en) 2017-06-14 2017-06-14 State estimation device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/021964 WO2018229898A1 (ja) 2017-06-14 2017-06-14 状態推定装置

Publications (1)

Publication Number Publication Date
WO2018229898A1 true WO2018229898A1 (ja) 2018-12-20

Family

ID=64660249

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/021964 WO2018229898A1 (ja) 2017-06-14 2017-06-14 状態推定装置

Country Status (4)

Country Link
US (1) US20200132775A1 (ja)
EP (1) EP3640822A4 (ja)
JP (1) JP6768155B2 (ja)
WO (1) WO2018229898A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021014499A1 (ja) * 2019-07-19 2021-01-28 三菱電機株式会社 パラメータ同定装置、パラメータ同定方法およびコンピュータプログラム
CN112847334A (zh) * 2020-12-16 2021-05-28 北京无线电测量研究所 一种基于视觉伺服的机械臂目标跟踪方法
CN115210609A (zh) * 2020-02-21 2022-10-18 株式会社东京测振 推定装置、振动感测器系统、由推定装置执行的方法及程序
US20220353277A1 (en) * 2021-04-19 2022-11-03 Daegu Gyeongbuk Institute Of Science And Technology Apparatus and method for estimating system state

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110727196B (zh) * 2019-09-26 2021-09-17 南京航空航天大学 基于鲁棒滤波器的正线性网络控制系统故障检测方法
US20220017103A1 (en) * 2020-07-17 2022-01-20 Toyota Research Institute, Inc. Systems and methods for automatically generating solver code for nonlinear model predictive control solvers
CN112115419A (zh) * 2020-09-14 2020-12-22 深圳大学 系统状态估计方法、系统状态估计装置
CN112491393B (zh) * 2020-11-30 2024-05-10 长沙师范学院 一种基于观测噪声协方差矩阵未知的线性滤波方法
CN113343453B (zh) * 2021-05-28 2023-02-10 华南理工大学 一种基于小步长离散化的电力电子级联变换器的建模方法
CN113626983B (zh) * 2021-07-06 2022-09-13 南京理工大学 基于状态方程递推预测高炮射弹脱靶量的方法
CN113472318B (zh) * 2021-07-14 2024-02-06 青岛杰瑞自动化有限公司 一种顾及观测模型误差的分级自适应滤波方法及系统
CN114994601A (zh) * 2022-06-02 2022-09-02 合肥联睿微电子科技有限公司 基于距离测量的广义卡尔曼滤波定位方法及系统
CN117521018B (zh) * 2024-01-08 2024-03-26 鹏城实验室 基于扩展观测的融合估计方法、装置、设备及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005248946A (ja) 2004-03-02 2005-09-15 General Electric Co <Ge> ガスタービンエンジンのためのモデルベース制御システム及び方法
JP2013036987A (ja) * 2011-07-08 2013-02-21 Canon Inc 情報処理装置及び情報処理方法
JP2013058105A (ja) * 2011-09-08 2013-03-28 Kyushu Institute Of Technology 物体運動推定装置、物体運動推定方法及びプログラム

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005248946A (ja) 2004-03-02 2005-09-15 General Electric Co <Ge> ガスタービンエンジンのためのモデルベース制御システム及び方法
JP2013036987A (ja) * 2011-07-08 2013-02-21 Canon Inc 情報処理装置及び情報処理方法
JP2013058105A (ja) * 2011-09-08 2013-03-28 Kyushu Institute Of Technology 物体運動推定装置、物体運動推定方法及びプログラム

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
FROGERAIS, PAUL ET AL.: "Various Ways to Compute the Continuous-Discrete Extended Kalman Filter", IEEE TRANSACTIONS ON AUTOMATIC CONTROL, vol. 57, no. 4, April 2012 (2012-04-01), pages 1000 - 1004, XP011440599, ISSN: 1558-2523 *
See also references of EP3640822A4
TAKASHI YAHAGI: "Kalman Filter and Adaptive Signal Processing", December 2005, CORONA PUBLISHING, pages: 48 - 51

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021014499A1 (ja) * 2019-07-19 2021-01-28 三菱電機株式会社 パラメータ同定装置、パラメータ同定方法およびコンピュータプログラム
TWI736358B (zh) * 2019-07-19 2021-08-11 日商三菱電機股份有限公司 參數鑑別裝置、參數鑑別方法及電腦程式
JPWO2021014499A1 (ja) * 2019-07-19 2021-12-09 三菱電機株式会社 パラメータ同定装置、パラメータ同定方法およびコンピュータプログラム
CN114127643A (zh) * 2019-07-19 2022-03-01 三菱电机株式会社 参数同定装置、参数同定方法及计算机程序
JP7066065B2 (ja) 2019-07-19 2022-05-12 三菱電機株式会社 パラメータ同定装置、パラメータ同定方法およびコンピュータプログラム
CN114127643B (zh) * 2019-07-19 2024-03-29 三菱电机株式会社 参数同定装置、参数同定方法及存储介质
CN115210609A (zh) * 2020-02-21 2022-10-18 株式会社东京测振 推定装置、振动感测器系统、由推定装置执行的方法及程序
CN112847334A (zh) * 2020-12-16 2021-05-28 北京无线电测量研究所 一种基于视觉伺服的机械臂目标跟踪方法
CN112847334B (zh) * 2020-12-16 2022-09-23 北京无线电测量研究所 一种基于视觉伺服的机械臂目标跟踪方法
US20220353277A1 (en) * 2021-04-19 2022-11-03 Daegu Gyeongbuk Institute Of Science And Technology Apparatus and method for estimating system state

Also Published As

Publication number Publication date
JPWO2018229898A1 (ja) 2019-11-07
US20200132775A1 (en) 2020-04-30
EP3640822A1 (en) 2020-04-22
EP3640822A4 (en) 2020-07-08
JP6768155B2 (ja) 2020-10-14

Similar Documents

Publication Publication Date Title
WO2018229898A1 (ja) 状態推定装置
Ding et al. Combined parameter and output estimation of dual-rate systems using an auxiliary model
Jouffroy et al. Finite-time simultaneous parameter and state estimation using modulating functions
Kumar et al. Solution of Fokker-Planck equation by finite element and finite difference methods for nonlinear systems
US10620597B2 (en) Gray box model estimation for process controller
JP2002049409A (ja) プロセス制御システムにおける適応推定モデル
Sun et al. A computationally efficient norm optimal iterative learning control approach for LTV systems
JP2015521748A (ja) 入力信号を変換する方法
Luan et al. Higher order moment stability region for Markov jump systems based on cumulant generating function
JP2019207660A (ja) 異常検出装置、異常検出方法、異常検出プログラム及び記録媒体
Zhang et al. Identification of continuous-time nonlinear systems: The nonlinear difference equation with moving average noise (NDEMA) framework
KR20200084010A (ko) 목표 시스템에 대한 제어 시스템 생성
JP2014041547A (ja) 時系列データ解析装置、方法、及びプログラム
JP2015219714A (ja) 状態推定装置、方法、及びプログラム
Garcia et al. Parameter estimation in time-triggered and event-triggered model-based control of uncertain systems
JP4815391B2 (ja) モデルパラメータ推定演算装置及び方法、モデルパラメータ推定演算処理プログラム並びにそれを記録した記録媒体
Üstüntürk Output feedback stabilization of nonlinear dual-rate sampled-data systems via an approximate discrete-time model
Gräßle Adaptivity in model order reduction with proper orthogonal decomposition
JP2022070386A (ja) 学習方法、系列解析方法、学習装置、系列解析装置、及びプログラム
Wang Optimized equivalent linearization for random vibration
Abgrall et al. A hybrid finite element–finite volume method for conservation laws
Jaramillo et al. Impulsive observer design for a class of nonlinear Lipschitz systems with time-varying uncertainties
Larsson et al. Estimation of continuous-time stochastic system parameters
Sfaihi et al. On stability conditions of singularly perturbed nonlinear Lur’e discrete-time systems
Dutra Uncertainty estimation in equality-constrained MAP and maximum likelihood estimation with applications to system identification and state estimation

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: 17913138

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2019524625

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: 2017913138

Country of ref document: EP

Effective date: 20200114