知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024054650
(43)【公開日】2024-04-17
(54)【発明の名称】処理装置、処理方法、およびプログラム
(51)【国際特許分類】
G01L 1/00 20060101AFI20240410BHJP
G01M 17/10 20060101ALI20240410BHJP
【ＦＩ】
G01L1/00 M
G01M17/10
【審査請求】未請求
【請求項の数】9
【出願形態】ＯＬ
(21)【出願番号】P 2022161018
(22)【出願日】2022-10-05
(71)【出願人】
【識別番号】000006655
【氏名又は名称】日本製鉄株式会社
(74)【代理人】
【識別番号】100090273
【弁理士】
【氏名又は名称】國分 孝悦
(72)【発明者】
【氏名】木村 幸彦
(72)【発明者】
【氏名】下川 嘉之
(72)【発明者】
【氏名】岸野 瞬士
(72)【発明者】
【氏名】伊藤 星太
(72)【発明者】
【氏名】山崎 陽介
(72)【発明者】
【氏名】藤本 隆裕
(72)【発明者】
【氏名】山尾 仁志
(72)【発明者】
【氏名】近藤 修
(72)【発明者】
【氏名】中川 淳一
(72)【発明者】
【氏名】南 秀樹
(72)【発明者】
【氏名】五百旗頭 学
(57)【要約】
【課題】 物体の応力を推定するに際し、計算速度を向上させる。
【解決手段】 処理装置４００は、軸箱１７ａ、１７ｂの加速度データを取得する。処理装置４００は、軸箱１７ａ、１７ｂの加速度データに基づいて、台車枠１６の応力に対する影響因子を算出する。処理装置４００は、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式に、算出した影響因子を与えることにより、台車枠１６の応力を算出する。
【選択図】 図１
【解決手段】 処理装置４００は、軸箱１７ａ、１７ｂの加速度データを取得する。処理装置４００は、軸箱１７ａ、１７ｂの加速度データに基づいて、台車枠１６の応力に対する影響因子を算出する。処理装置４００は、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式に、算出した影響因子を与えることにより、台車枠１６の応力を算出する。
【選択図】 図１
【特許請求の範囲】
【請求項1】
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得手段と、
前記データ取得手段により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出手段と、
を有し、
前記応力算出手段は、前記データ取得手段により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理装置。
前記データ取得手段により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出手段と、
を有し、
前記応力算出手段は、前記データ取得手段により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理装置。
【請求項2】
前記複数の説明変数は、変換前の第２物理量と、当該変換前の第２物理量に乗算される係数と、に基づいて、変換後の前記第２物理量を算出する計算式における、前記係数、または、前記変換後の第２物理量を含み、
前記第２物理量は、位置、変位、および速度のうち、少なくとも１つを含む、請求項１に記載の処理装置。
前記第２物理量は、位置、変位、および速度のうち、少なくとも１つを含む、請求項１に記載の処理装置。
【請求項3】
前記計算式は、前記物体の変位前の物理座標系の位置の座標と、当該座標に乗算される変換行列と、に基づいて、前記物体の変位後の物理座標系の位置を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、請求項２に記載の処理装置。
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、請求項２に記載の処理装置。
【請求項4】
前記計算式は、物理座標系における前記物体の振動を記述する運動方程式を固有値解析することにより算出される固有ベクトルと、当該固有ベクトルに乗算されるモード座標系における前記物体の変位と、に基づいて、物理座標系における前記物体の変位を算出する計算式と、前記固有ベクトルと、モード座標系における前記物体の速度と、に基づいて、物理座標系における前記物体の速度を算出する計算式と、を含み、
前記複数の説明変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、請求項２に記載の処理装置。
前記複数の説明変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、請求項２に記載の処理装置。
【請求項5】
前記計算式は、前記物体のモード座標系における前記物体の変位と、当該座標に乗算される変換行列と、に基づいて、当該変位が生じている時刻の次の時刻での、前記物体のモード座標系における前記物体の変位を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、請求項２に記載の処理装置。
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、請求項２に記載の処理装置。
【請求項6】
前記応力算出手段は、モード座標系における前記物体の振動を表す運動方程式に基づいて構成される状態方程式と、観測変数と状態変数との関係を表す観測方程式と、前記第１物理量の測定データから算出される前記観測変数の測定値と、に基づいて、データ同化を行うフィルタを用いた演算を行うことにより、前記状態変数の推定値を決定し、
前記状態変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、請求項４または５に記載の処理装置。
前記状態変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、請求項４または５に記載の処理装置。
【請求項7】
前記物体は、鉄道車両の台車枠を含む、請求項１～５のいずれか１項に記載の処理装置。
【請求項8】
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得工程と、
前記データ取得工程により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出工程と、
を有し、
前記応力算出工程は、前記データ取得工程により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理方法。
前記データ取得工程により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出工程と、
を有し、
前記応力算出工程は、前記データ取得工程により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理方法。
【請求項9】
請求項１～５のいずれか１項に記載の処理装置の各手段としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、処理装置、処理方法、およびプログラムに関する。
【背景技術】
【0002】
鉄道台車の台車枠等、外力が作用する物体の応力を推定する技術が求められる。この種の技術として、特許文献１、２に記載の技術がある。
特許文献１には、以下の技術が開示されている。まず、変位前の位置座標と変位後の位置座標とに基づいてアフィン変換行列（変位前の位置座標に乗算される行列）を算出し、アフィン変換行列を特異値分解して特異値を対角成分として有する対角行列を算出する。そして、対角行列に基づいて主歪みを算出し、主歪みに基づいて応力を算出する。
特許文献１には、以下の技術が開示されている。まず、変位前の位置座標と変位後の位置座標とに基づいてアフィン変換行列（変位前の位置座標に乗算される行列）を算出し、アフィン変換行列を特異値分解して特異値を対角成分として有する対角行列を算出する。そして、対角行列に基づいて主歪みを算出し、主歪みに基づいて応力を算出する。
【0003】
また、特許文献２には、以下の技術が開示されている。まず、モード座標系における状態変数の時間変化を規定するための状態方程式と、観測変数および状態変数を関連付けるための観測方程式と、観測変数の測定値と、に基づいて、データ同化を行うフィルタを用いた演算を行うことにより、モード座標系における状態変数の推定値を算出する。そして、モード座標系における状態変数の推定値に基づいて、モード座標系における台車枠の変位分布を算出し、モード座標系における台車枠の変位分布を物理座標系の台車枠の変位分布に変換する。そして、物理座標系の台車枠の変位分布に基づいて、台車枠の歪み分布を算出し、台車枠の歪み分布に基づいて、台車枠の応力分布を算出する。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開２０２１－１９３３４６号公報
【特許文献2】特開２０２２－２８３７４号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかしながら、特許文献１、２に記載の技術では、物体の応力をリアルタイムで算出することが容易ではない。例えば、物体の応力は外力によって時々刻々と変化し得るところ、特許文献１に記載の技術では、応力の算出タイミングのそれぞれにおいて、行列の特異値分解を行う必要がある。また、特許文献２に記載の技術では、応力の算出タイミングのそれぞれにおいて、モード座標系における状態変数の推定値を導出したうえで、モード座標系における変位を物理座標系における変位に変換する必要がある。また、特許文献１、２に記載の技術では、台車枠の応力を推定する際に、台車枠の変位に基づいて、台車枠の歪みを算出し、台車枠の歪みに基づいて、台車枠の応力を算出する必要がある。以上のように、特許文献１、２に記載の技術では、物体の応力を高速に推定することができないという問題点がある。
【0006】
本発明は、以上のような問題点に鑑みてなされたものであり、物体の応力を推定するに際し、計算速度を向上させることを目的とする。
【課題を解決するための手段】
【0007】
本発明の処理装置は、物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得手段と、前記データ取得手段により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出手段と、を有し、前記応力算出手段は、前記データ取得手段により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する。
【0008】
本発明の処理方法は、物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得工程と、前記データ取得工程により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出工程と、を有し、前記応力算出工程は、前記データ取得工程により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する。
【0009】
本発明のプログラムは、前記処理装置の各手段としてコンピュータを機能させる。
【発明の効果】
【0010】
本発明によれば、物体の応力を推定するに際し、計算速度を向上させることができる。
【図面の簡単な説明】
【0011】
【図1】鉄道車両の概略の一例を示す図である。
【図2】台車枠およびその周辺の部品の構成の一例を示す図である。
【図3】連結器の一例をモデル化して示す図である。
【図4】処理装置の機能的な構成の一例を示す図である。
【図5】線形回帰式の回帰係数を算出する方法の一例を説明するフローチャートである。
【図6】台車枠の応力を算出する方法の一例を説明するフローチャートである。
【図7】台車枠の最大主応力と時間との関係の一例を示す図である。
【図8】第１実施形態の計算例を示し、台車枠の最大主応力の真値と推定値との関係を示す図である。
【図9】第２実施形態の計算例を示し、台車枠の最大主応力の真値と推定値との関係の第１の例を示す図である。
【図10】第２実施形態の計算例を示し、台車枠の最大主応力の真値と推定値との関係の第２の例を示す図である。
【図11】第２の実施形態の変形例の計算例を示し、台車枠の応力測定値と時間との関係の一例を示す図である。
【図12】第２の実施形態の変形例の計算例を示し、台車枠の応力の真値と推定値との関係の第１の例を示す図である。
【図13】第２の実施形態の変形例の計算例を示し、バネ帽１位の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図14】第２の実施形態の変形例の計算例を示し、バネ帽２位の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図15】第２の実施形態の変形例の計算例を示し、バネ帽３位の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図16】第２の実施形態の変形例の計算例を示し、バネ帽４位の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図17】第２の実施形態の変形例の計算例を示し、モータ座１軸の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図18】第２の実施形態の変形例の計算例を示し、モータ座２軸の上下方向の変位の測定値と計算値との関係の一例を示す図である。
【図19】第２の実施形態の変形例の計算例を示し、台車枠の応力の真値と推定値との関係の第２の例を示す図である。
【発明を実施するための形態】
【0012】
以下、図面を参照しながら、本発明の実施形態を説明する。
なお、長さ、位置、大きさ、間隔等、比較対象が同じであることは、厳密に同じである場合の他、発明の主旨を逸脱しない範囲で異なるもの（例えば、設計時に定められる公差の範囲内で異なるもの）も含むものとする。
なお、長さ、位置、大きさ、間隔等、比較対象が同じであることは、厳密に同じである場合の他、発明の主旨を逸脱しない範囲で異なるもの（例えば、設計時に定められる公差の範囲内で異なるもの）も含むものとする。
【0013】
以下の各実施形態の処理装置は、物体に生じる応力を算出することを含む処理を行う。外力により応力が生じる物体であれば処理装置における応力の算出対象の物体は限定されないが、以下の各実施形態では、鉄道車両の台車枠に生じる応力を処理装置が算出する場合を例示する。また、以下の各実施形態では、鉄道車両の台車枠の運動方程式を用いて、台車枠の変位を算出する場合を例示する。そこで、実施形態の説明の前に、鉄道車両の構成の一例と、台車枠の物理座標系における運動方程式の一例と、を概説する。
【0014】
（鉄道車両の概略構成）
図１は、鉄道車両の概略の一例を示す図である。図２は、台車枠およびその周辺の部品の構成の一例を示す図である。なお、図１および図２において、鉄道車両は、ｘ軸の正の方向に進むものとする（ｘ軸は、鉄道車両の走行方向に沿う軸である）。また、ｚ軸は、軌道２０（地面）に対し垂直方向（鉄道車両の高さ方向）であるものとする。ｙ軸は、鉄道車両の走行方向に対して垂直な水平方向（鉄道車両の走行方向と高さ方向との双方に垂直な方向）であるものとする。また、ｘ－ｙ－ｚ座標は、物理座標系における座標であるものとする。なお、鉄道車両は、営業車両であっても、試験用車両であっても、検査用車両であってもよい。また、各図において、○（白丸）の中に●（黒丸）が付されている記号は、紙面の奥側から手前側に向かう矢印線を示す。
図１は、鉄道車両の概略の一例を示す図である。図２は、台車枠およびその周辺の部品の構成の一例を示す図である。なお、図１および図２において、鉄道車両は、ｘ軸の正の方向に進むものとする（ｘ軸は、鉄道車両の走行方向に沿う軸である）。また、ｚ軸は、軌道２０（地面）に対し垂直方向（鉄道車両の高さ方向）であるものとする。ｙ軸は、鉄道車両の走行方向に対して垂直な水平方向（鉄道車両の走行方向と高さ方向との双方に垂直な方向）であるものとする。また、ｘ－ｙ－ｚ座標は、物理座標系における座標であるものとする。なお、鉄道車両は、営業車両であっても、試験用車両であっても、検査用車両であってもよい。また、各図において、○（白丸）の中に●（黒丸）が付されている記号は、紙面の奥側から手前側に向かう矢印線を示す。
【0015】
図１および図２に示す例では、鉄道車両は、車体１１と、台車１２ａ、１２ｂと、輪軸１３ａ～１３ｄと、を有する。このように図１および図２では、１つの車体１１に、２つの台車１２ａ、１２ｂと、４組の輪軸１３ａ～１３ｄと、が備わる鉄道車両を例示する。輪軸１３ａ～１３ｄは、車軸１５ａ～１５ｄと、その両端に設けられた車輪１４ａ～１４ｄと、を有する。図１および図２では、台車１２ａ、１２ｂが、ボルスタレス台車である場合を例示する。
【0016】
図１では、表記の都合上、輪軸１３ａ～１３ｄの一方の車輪１４ａ～１４ｄのみを示すが、輪軸１３ａ～１３ｄの他方にも車輪が配置されている（図１に示す例では、車輪は合計８つある）。
図２において、輪軸１３ａ、１３ｂのｙ軸に沿う方向の両側には、軸箱１７ａ、１７ｂが配置される。軸箱１７ａ、１７ｂは、モノリンク１８ａ、１８ｂを介して台車枠１６と接続される。また、軸箱１７ａ、１７ｂは、軸バネ１９ａ、１９ｂを介して台車枠１６と接続される。なお、鉄道車両は、図１および図２に示す構成要素以外の構成要素を有する。表記および説明の都合上、図１および図２では、当該構成要素の図示を省略する。例えば、鉄道車両が軸ダンパを有する場合、軸箱１７ａ、１７ｂは、軸ダンパを介して台車枠１６と接続されてもよい。
図２において、輪軸１３ａ、１３ｂのｙ軸に沿う方向の両側には、軸箱１７ａ、１７ｂが配置される。軸箱１７ａ、１７ｂは、モノリンク１８ａ、１８ｂを介して台車枠１６と接続される。また、軸箱１７ａ、１７ｂは、軸バネ１９ａ、１９ｂを介して台車枠１６と接続される。なお、鉄道車両は、図１および図２に示す構成要素以外の構成要素を有する。表記および説明の都合上、図１および図２では、当該構成要素の図示を省略する。例えば、鉄道車両が軸ダンパを有する場合、軸箱１７ａ、１７ｂは、軸ダンパを介して台車枠１６と接続されてもよい。
【0017】
図２では、１つの台車１２ａに１つの台車枠１６が配置される場合を例示する。軸箱１７ａ、１７ｂ、モノリンク１８ａ、１８ｂ、および軸バネ１９ａ、１９ｂは、１つの車輪に対して１つずつ配置される。前述したように１つの台車１２ａには、４つの車輪が配置される。したがって、１つの台車１２ａには、軸箱、モノリンク、および軸バネが、それぞれ４つずつ配置される。
【0018】
図２では、台車１２ａの構成（台車枠１６、軸箱１７ａ、１７ｂ、モノリンク１８ａ、１８ｂ、および軸バネ１９ａ、１９ｂ）のみを示す。台車１２ｂの構成（台車枠、軸箱、モノリンク、および軸バネ）も図２に示す構成と同じ構成で実現される。なお、鉄道車両自体は公知の技術で実現できる。したがって、ここでは、鉄道車両の概要のみを説明し、その詳細な説明を省略する。また、鉄道車両は、図１および図２に示すものに限定されない。
以下の説明では、台車枠１６と、軸箱１７ａ、１７ｂと、を連結する部品を、必要に応じて連結器と総称する。
以下の説明では、台車枠１６と、軸箱１７ａ、１７ｂと、を連結する部品を、必要に応じて連結器と総称する。
【0019】
図１および図２では、軸箱１７ａ、１７ｂのそれぞれに加速度センサ２１ａ、２１ｂが取り付けられる場合を例示する。また、図１および図２では、台車枠１６にも加速度センサ２２ａ、２２ｂが取り付けられる場合を例示する。また、図１および図２では、加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂが、３次元加速度センサである場合を例示する。加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂで測定される加速度のデータから、加速度のｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分が得られる。また、加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂで測定された加速度のデータから、加速度のｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分について、時間に関する２階積分を行うことにより、加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂの取り付け位置での変位が算出される。加速度センサ２１ａ、２１ｂは、軸箱１７ａ、１７ｂの所定の位置に取り付けられ、加速度センサ２２ａ、２２ｂは、台車枠１６の所定の位置に取り付けられる。図２では、軸箱１７ａ、１７ｂに取り付けられる加速度センサ２１ａ、２１ｂの数が１であり、台車枠１６に取り付けられる加速度センサ２２ａ、２２ｂの数が２である場合を例示する。しかしながら、軸箱１７ａ、１７ｂに取り付けられる加速度センサ２１ａ、２１ｂの数、および、台車枠１６に取り付けられる加速度センサの数は、１以上であれば、幾つであってもよい。
【0020】
連結器は、台車枠１６および軸箱１７ａ、１７ｂに接続される。したがって、軸箱１７ａ、１７ｂの振動と連結器の振動とは連動する。これらの振動は台車枠１６に伝搬する。台車枠１６に作用する外力は、連結器の振動における粘性減衰力と剛性力との和で表される。粘性減衰力は、粘性減衰係数と速度との積で表される。剛性力は、剛性と変位との積で表される。
【0021】
台車枠１６に作用する外力が得られれば、台車枠１６の振動（運動）を表す運動方程式を用いて、台車枠１６の変位（台車枠１６の各部の位置の変化量）を算出することができる。また、台車枠１６の変位を空間微分（偏微分）することにより、台車枠１６の各種の物理量（例えば、歪み、応力、および表面力）を算出することができる。
【0022】
（台車枠１６の物理座標系における運動方程式）
台車枠１６の振動を表す物理座標系における運動方程式は、例えば以下の（１）式で表される。なお、後述する第１実施形態では、（１）式で表される物理座標系における運動方程式を用いる。一方、後述する第２実施形態では、物理座標系における運動方程式に対応するモード座標系における運動方程式を用いる（（５１）式を参照）。
台車枠１６の振動を表す物理座標系における運動方程式は、例えば以下の（１）式で表される。なお、後述する第１実施形態では、（１）式で表される物理座標系における運動方程式を用いる。一方、後述する第２実施形態では、物理座標系における運動方程式に対応するモード座標系における運動方程式を用いる（（５１）式を参照）。
【0023】
【数1】
【0024】
ここで、［Ｍ］（∈Ｒ3l×3l）は、台車枠１６の質量行列である。［Ｃ］（∈Ｒ3l×3l）は、台車枠１６の粘性行列（減衰行列ともいう）である。［Ｋ］（∈Ｒ3l×3l）は、台車枠１６の剛性行列である。
｛ｕ｝（∈Ｒ3l）は、台車枠１６の変位ベクトルである。｛ｆ｝（∈Ｒ3l）は、台車枠１６の外力ベクトルである。
｛ｕ｝（∈Ｒ3l）は、台車枠１６の変位ベクトルである。｛ｆ｝（∈Ｒ3l）は、台車枠１６の外力ベクトルである。
【0025】
以下の各実施形態では、物理座標系において、台車枠１６の振動には、ｘ軸方向、ｙ軸方向、およびｚ軸方向の３つの成分がある場合を例示する。したがって、｛ｕ｝を構成する各節点の変位および｛ｆ｝を構成する各節点の外力は、３つの自由度を有する。なお、（１）式において、・は、ｄ／ｄｔ（時間の１階微分）を表し、・・は、ｄ2／ｄｔ2（時間の２階微分）を表す（このことは、以降の式でも同じである）。
【0026】
ｌは、数値解析における変位の近似解の自由度に対応する。以下の各実施形態では、数値解析として、有限要素法（ＦＥＭ；Finite Element Method）を用いる場合を例示する。したがって、ｌは、例えば、有限要素法において定められる要素の節点の数である（以下の説明では要素（小領域）を必要に応じてメッシュと称する）。この場合、台車枠１６の質量行列［Ｍ］の成分、台車枠１６の粘性行列［Ｃ］の成分、台車枠１６の剛性行列［Ｋ］の成分には、それぞれ、各成分に対応する密度から算出される値、各成分に対応する粘性減衰係数から算出される値、各成分に対応する剛性から算出される値が与えられる。密度、粘性減衰係数、および剛性は、位置によらずに同じ値としてもよいし、異なる値としてもよい。台車枠１６の質量行列［Ｍ］の成分、粘性行列［Ｃ］の成分、および剛性行列［Ｋ］の成分は、例えば、有限要素法による数値解析を行う公知のソルバーにおいて、有限要素法のメッシュと、台車枠１６全体の密度、粘性減衰係数、および剛性と、を用いて算出される。なお、数値解析の手法は、有限要素法に限定されず、有限要素法以外のその他の手法（例えば、有限差分法（ＦＤＭ；Finite Difference Method））であってもよい。
【0027】
（１）式の左辺第１項は、台車枠１６に作用する重力を表す慣性項である。（１）式の左辺第２項は、台車枠１６に作用する粘性力を表す減衰項である。（１）式の左辺第３項は、台車枠１６に作用する剛性力を表す剛性項である。
【0028】
（１）式の右辺の外力ベクトル｛ｆ｝は、連結器に作用する外力を算出することによって与えられる。
そこで、図３を参照しながら、連結器に作用する外力の算出方法の一例を説明する。図３は、連結器の一例をモデル化して示す図である。図３（ａ）は、台車枠１６と軸箱１７ａとに接続される連結器をモデル化した図の一例を示す。図３（ｂ）は、台車枠１６と軸箱１７ｂとに接続される連結器をモデル化した図の一例を示す。台車枠１６とその他の軸箱とに接続される連結器をモデル化した図も、図３（ａ）および図３（ｂ）と同じようにして表されるので、ここでは、その詳細な説明を省略する。台車枠１６と連結器との接続箇所は、台車枠１６と連結器とが相互に接触する領域（の全体）としても、台車枠１６と連結器とが接触する領域の代表点（例えば、重心の位置）としてもよい。図３に示す例では、説明を簡単にするため、台車枠１６と連結器との接続箇所が点であるものとする。以下の説明では、台車枠１６と連結器との接続箇所を、必要に応じて、着力点と称する。図３に示す例では、バネおよびダンパを並列に接続したモデルで、モノリンク１８ａおよび軸バネ１９ａを表す。モデル化したモノリンク１８ａと台車枠１６とは着力点３１ａで接続される。また、モデル化した軸バネ１９ａと台車枠１６とは着力点３２ａで接続される。また、軸箱１７ｂについても同様に、モデル化したモノリンク１８ｂと台車枠１６とは着力点３１ｂで接続され、モデル化した軸バネ１９ｂと台車枠１６とは着力点３２ｂで接続される。
連結器の振動を表す物理座標系における運動方程式は、以下の（２）式で表される。
そこで、図３を参照しながら、連結器に作用する外力の算出方法の一例を説明する。図３は、連結器の一例をモデル化して示す図である。図３（ａ）は、台車枠１６と軸箱１７ａとに接続される連結器をモデル化した図の一例を示す。図３（ｂ）は、台車枠１６と軸箱１７ｂとに接続される連結器をモデル化した図の一例を示す。台車枠１６とその他の軸箱とに接続される連結器をモデル化した図も、図３（ａ）および図３（ｂ）と同じようにして表されるので、ここでは、その詳細な説明を省略する。台車枠１６と連結器との接続箇所は、台車枠１６と連結器とが相互に接触する領域（の全体）としても、台車枠１６と連結器とが接触する領域の代表点（例えば、重心の位置）としてもよい。図３に示す例では、説明を簡単にするため、台車枠１６と連結器との接続箇所が点であるものとする。以下の説明では、台車枠１６と連結器との接続箇所を、必要に応じて、着力点と称する。図３に示す例では、バネおよびダンパを並列に接続したモデルで、モノリンク１８ａおよび軸バネ１９ａを表す。モデル化したモノリンク１８ａと台車枠１６とは着力点３１ａで接続される。また、モデル化した軸バネ１９ａと台車枠１６とは着力点３２ａで接続される。また、軸箱１７ｂについても同様に、モデル化したモノリンク１８ｂと台車枠１６とは着力点３１ｂで接続され、モデル化した軸バネ１９ｂと台車枠１６とは着力点３２ｂで接続される。
連結器の振動を表す物理座標系における運動方程式は、以下の（２）式で表される。
【0029】
【数2】
【0030】
ここで、［Ｃbc］は、連結器の粘性行列である。［Ｋbc］は、連結器の剛性行列である。｛ｕ0｝は、連結器の軸箱１７ａとの接続箇所における変位と、連結器の着力点における変位と、で構成される連結器の変位ベクトルである。｛ｆ0｝は、連結器の軸箱１７ａとの接続箇所における外力と、連結器の着力点における外力と、で構成される連結器の外力ベクトルである。連結器の粘性行列［Ｃbc］の成分、連結器の剛性行列［Ｋbc］の成分には、それぞれ、各成分に対応する粘性減衰係数から算出される値、各成分に対応する剛性から算出される値が与えられる。連結器の粘性行列［Ｃ］および剛性行列［Ｋ］の成分は、例えば、連結器の軸箱１７ａとの接続箇所の位置と、連結器の着力点の位置と、連結器の密度と、連結器の粘性減衰係数と、連結器の剛性と、を用いて算出される。
【0031】
（２）式において、着力点におけるｘ軸方向、ｙ軸方向、およびｚ軸方向の変位を０（零）とする。また、連結器の軸箱１７ａとの接続箇所における変位を軸箱１７ａの変位とする。以上のようにすることによって、連結器の外力ベクトル｛ｆ0｝を構成する外力として、連結器の着力点に作用する外力を算出することができる。つまり、モノリンク１８ａの着力点３１ａおよび軸バネ１９ａの着力点３２ａに作用する外力は、モノリンク１８ａおよび軸バネ１９ａのそれぞれについて構成した（２）式を用いて算出される。また、軸箱１７ｂについても同様に、軸箱１７ｂの変位から、モノリンク１８ｂの着力点３１ｂおよび軸バネ１９ｂの着力点３２ｂに作用する外力が算出される。そして、台車枠１６の着力点３１ａ、３１ｂに作用する外力、着力点３２ａ、３２ｂに作用する外力は、それぞれ、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力の反作用力、軸バネ１９ａ、１９ａの着力点３２ａ、３２ｂに作用する外力の反作用力として算出される。
【0032】
（１）式の外力ベクトル｛ｆ｝の成分のうち、（台車枠１６の）着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力の成分には、前述の方法で算出した値を与え、その他の成分には、０（零）を与えて、台車枠１６の変位ベクトル｛ｕ｝を算出することにより、台車枠１６の変位を算出することができる。
なお、台車枠１６の運動方程式自体は公知の方程式でよく、前述した運動方程式に限定されない。また、台車枠１６の運動方程式から台車枠１６の変位を算出する手法自体は、公知の技術で実現され、前述した手法に限定されない。
なお、台車枠１６の運動方程式自体は公知の方程式でよく、前述した運動方程式に限定されない。また、台車枠１６の運動方程式から台車枠１６の変位を算出する手法自体は、公知の技術で実現され、前述した手法に限定されない。
【0033】
（第１実施形態）
次に、第１実施形態を説明する。
＜処理装置４００、処理方法＞
図４は、処理装置４００の機能的な構成の一例を示す図である。処理装置４００のハードウェアは、例えば、ＣＰＵ（Central Processing Unit）等、一またはそれ以上の数のハードウェアプロセッサを有する。また、処理装置４００のハードウェアは、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等、一またはそれ以上の数のメモリを有する。処理装置４００は、メモリに格納される一またはそれ以上の数のプログラムを一またはそれ以上の数のハードウェアプロセッサにより実行することで各種の演算を実行する。さらに、処理装置４００のハードウェアは、入力装置、および出力装置を有する。また、処理装置４００は、ＡＳＩＣ（Application Specific Integrated Circuit）等の専用のハードウェアにより実現されてもよい。また、図１に示すように本実施形態では、処理装置４００が鉄道車両の車体１１内に設置される場合を例示する。しかしながら、処理装置４００は鉄道車両の外部に設置されてもよい。
次に、第１実施形態を説明する。
＜処理装置４００、処理方法＞
図４は、処理装置４００の機能的な構成の一例を示す図である。処理装置４００のハードウェアは、例えば、ＣＰＵ（Central Processing Unit）等、一またはそれ以上の数のハードウェアプロセッサを有する。また、処理装置４００のハードウェアは、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等、一またはそれ以上の数のメモリを有する。処理装置４００は、メモリに格納される一またはそれ以上の数のプログラムを一またはそれ以上の数のハードウェアプロセッサにより実行することで各種の演算を実行する。さらに、処理装置４００のハードウェアは、入力装置、および出力装置を有する。また、処理装置４００は、ＡＳＩＣ（Application Specific Integrated Circuit）等の専用のハードウェアにより実現されてもよい。また、図１に示すように本実施形態では、処理装置４００が鉄道車両の車体１１内に設置される場合を例示する。しかしながら、処理装置４００は鉄道車両の外部に設置されてもよい。
【0034】
本実施形態では、処理装置４００が、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式を用いて、台車枠１６の応力を算出する場合を例示する。したがって、本実施形態では、線形回帰式の回帰係数を算出する学習と、学習により算出された回帰係数が与えられた線形回帰式を用いて台車枠１６の応力を算出する推定と、が行われる。図５は、線形回帰式の回帰係数を算出する方法の一例を説明するフローチャートである。図６は、台車枠１６の応力を算出する方法の一例を説明するフローチャートである。図５および図６のフローチャートは、例えば、処理装置４００のハードウェアプロセッサが、処理装置４００のメモリに格納される一またはそれ以上の数のプログラムを実行することで実現される。なお、説明変数は独立変数等とも称され、目的変数は従属変数等とも称される。また、本実施形態の線形回帰式は、複数の説明変数を含むので重回帰式である。
【0035】
<<データ取得部４１０、ステップＳ５０１、Ｓ６０１>>
データ取得部４１０は、処理装置４００における計算のために予め処理装置４００が取得しておく必要がある各種の情報を取得する。以下に本実施形態のデータ取得部４１０が取得する情報の一例を説明する。
データ取得部４１０は、処理装置４００における計算のために予め処理装置４００が取得しておく必要がある各種の情報を取得する。以下に本実施形態のデータ取得部４１０が取得する情報の一例を説明する。
【0036】
データ取得部４１０は、加速度センサ２１ａ、２１ｂで測定された軸箱１７ａ、１７ｂの加速度のデータを取得する。以下の説明では、加速度のデータを必要に応じて加速度データと称する。また、データ取得部４１０は、加速度センサ２２ａ、２２ｂで測定された台車枠１６の加速度データを取得する。
【0037】
本実施形態では、処理装置４００（データ取得部４１０）と、加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂと、が有線または無線により通信可能に接続される場合を例示する。また、本実施形態では、データ取得部４１０が、加速度センサ２１ａ～２１ｂ、２２ａ～２２ｂから、軸箱１７ａ、１７ｂの加速度データ、台車枠１６の加速度データをそれぞれ取得する場合を例示する。しかしながら、必ずしもこのようにする必要はない。例えば、処理装置４００（データ取得部４１０）は、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データを、加速度センサ２１ａ、２１ｂ、２２ａ、２２ｂとは異なる外部装置を経由して取得してもよい。本実施形態では、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データが、物体の振動に応じて変化する第１物理量の測定データである場合を例示する。
【0038】
データ取得部４１０は、解析モデルを取得する。前述したように本実施形態では、数値解析として、有限要素法を用いる場合を例示する。したがって本実施形態では、データ取得部４１０は、例えば、台車枠１６および連結器を含む領域の形状を複数のメッシュで表した場合の各メッシュの情報を含む情報を、解析モデルの情報として取得する。解析モデルの作成自体は、公知の手法で実現されるので、ここでは解析モデルの詳細な説明を省略する。また、データ取得部４１０は、各メッシュのうち、後述する部位ｉ、ｊの位置に対応するメッシュの情報を、解析モデルの情報として取得する。
【0039】
データ取得部４１０は、処理装置４００における計算の際に定数として扱われる情報を取得する。データ取得部４１０は、例えば、台車枠１６および連結器のそれぞれについて、密度、粘性減衰係数、および剛性を取得する。また、データ取得部４１０は、例えば、台車枠１６のヤング率およびポアソン比を取得する。これらの情報は、例えば、解析モデル（各メッシュ）に与えられる。
【0040】
データ取得部４１０が、解析モデルの情報と、処理装置４００における計算の際に定数として扱われる情報と、を取得する方法は限定されない。データ取得部４１０は、例えば、処理装置４００の入力装置に対するオペレータの情報入力操作に基づいてこれらの情報を入力してもよいし、外部装置からこれらの情報を受信してもよいし、処理装置４００に接続された記憶媒体からこれらの情報を読み出してもよい。
【0041】
なお、ステップＳ５０１、Ｓ６０１においては、データ取得部４１０が取得するデータの時刻が異なる。ステップＳ５０１においては、学習時刻ｔ１におけるデータが取得される。一方、ステップＳ６０１においては、学習時刻ｔ１よりも後の時刻である推定時刻ｔ２におけるデータが取得される。ただし、学習時刻ｔ１および推定時刻ｔ２で共通するデータについては、ステップＳ５０１で取得されたデータを記憶しておき、ステップＳ６０１で読み出してもよい。
【0042】
<<回帰係数算出部４２０、回帰係数記憶部４３０、ステップＳ５０２、Ｓ５０３>>
回帰係数算出部４２０は、データ取得部４１０により取得されたデータを用いて、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式の回帰係数を算出する。本実施形態では、台車枠１６の応力（目的変数）が、台車枠１６の最大主応力である場合を例示する。
回帰係数算出部４２０は、データ取得部４１０により取得されたデータを用いて、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式の回帰係数を算出する。本実施形態では、台車枠１６の応力（目的変数）が、台車枠１６の最大主応力である場合を例示する。
【0043】
また、本実施形態では、複数の説明変数が、台車枠１６の部位ｉの座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）をアフィン変換する際に当該座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）に乗算されるアフィン変換行列の成分のうち、線形変換を行うための成分である場合を例示する。台車枠１６の部位ｉの座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）をアフィン変換することにより、台車枠１６の部位ｉの座標が（ｘ2
(i)，ｙ2
(i)，ｚ2
(i)）に変化することが表される。
【0044】
ここで、線形変換を行うための成分は、平行移動を行うための成分を除く成分である。線形変換を行うための成分は、拡縮（拡大および縮小）、回転、およびせん断のうちの少なくとも１つを行うための成分である。なお、ｘ1
(i)，ｙ1
(i)，ｚ1
(i)、ｘ2
(i)，ｙ2
(i)，ｚ2
(i)の、１、２は、それぞれ、変換前の座標、変換後の座標であることを示す。また、各変数において上付きの（ｉ）は、部位ｉの値であることを表す記号である。台車枠１６の部位ｉは、例えば、メッシュの中から選ばれる。
【0045】
以下に本実施形態における線形回帰式の一例を説明する。
アフィン変換係数をａ1～ａ12とする。ｎ1個の部位ｉ（ｉ＝１，２，・・・，ｎ1）の座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）は、以下の（３）式のように、座標（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）にアフィン変換される。ここで、ｎ1≧４（ｎ1が４以上の整数であること）が必要条件である。ｎ1＝４であれば、アフィン変換係数ａ1～ａ12を算出することができる。ｎ1の数を多くすれば、アフィン変換係数ａ1～ａ12の計算精度は高まるが計算負荷も高まる。このような観点からｎ1は予め設定される。また、部位ｉは、台車枠１６の応力（本実施形態では最大主応力）の算出対象の位置に対応する部位ｊに応じて予め定められる。例えば、部位ｉは、ｎ1個の部位ｉのうちの任意の４つの部位ｉを結ぶ線を境界線とする領域の内側に、台車枠１６の応力の算出対象の位置に対応する部位ｊが位置するように定められてもよい。また、部位ｉは、例えば、台車枠１６の応力の算出対象の位置に対応する部位ｊの近傍に（すなわち、部位ｉ、ｊの距離が所定距離以下となるように）定められてもよい。
アフィン変換係数をａ1～ａ12とする。ｎ1個の部位ｉ（ｉ＝１，２，・・・，ｎ1）の座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）は、以下の（３）式のように、座標（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）にアフィン変換される。ここで、ｎ1≧４（ｎ1が４以上の整数であること）が必要条件である。ｎ1＝４であれば、アフィン変換係数ａ1～ａ12を算出することができる。ｎ1の数を多くすれば、アフィン変換係数ａ1～ａ12の計算精度は高まるが計算負荷も高まる。このような観点からｎ1は予め設定される。また、部位ｉは、台車枠１６の応力（本実施形態では最大主応力）の算出対象の位置に対応する部位ｊに応じて予め定められる。例えば、部位ｉは、ｎ1個の部位ｉのうちの任意の４つの部位ｉを結ぶ線を境界線とする領域の内側に、台車枠１６の応力の算出対象の位置に対応する部位ｊが位置するように定められてもよい。また、部位ｉは、例えば、台車枠１６の応力の算出対象の位置に対応する部位ｊの近傍に（すなわち、部位ｉ、ｊの距離が所定距離以下となるように）定められてもよい。
【0046】
【数3】
【0047】
（３）式は、部位ｉにおいて弾性変形等が生じると、部位ｉの座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）は、座標（ｘ2
(i)，ｙ2
(i)，ｚ2
(i)）に変化することを示す。
ここで、係数行列Ａ、変換前座標行列Ｃ1、変換後座標行列Ｃ2を、それぞれ、以下の（４）式、（５）式、（６）式のように定義する。
ここで、係数行列Ａ、変換前座標行列Ｃ1、変換後座標行列Ｃ2を、それぞれ、以下の（４）式、（５）式、（６）式のように定義する。
【0048】
【数4】
【0049】
係数行列Ａは、重回帰分析における教師あり学習を行うことにより算出される。ここでは、最小二乗法により係数行列Ａを算出する場合を例示する。この場合、係数行列Ａは、例えば、以下の（７）式により算出される。なお、（７）式は、（７）式の［］内の値が最小になるときの係数行列Ａを算出することを示す数式である。また、（７）式のｃ1,i、ｃ2,iは、それぞれ、以下の（８）式、（９）式で表される。
【0050】
【数5】
【0051】
ここで、（７）式の［］内の値が最小になるときの係数行列Ａを係数行列Ａ＾と表記する。また、係数行列Ａ＾の成分であるアフィン係数をａ1＾、ａ2＾、・・・、ａ12＾と表記する（以下の（１０）式を参照）。なお、Ａ＾、ａ1＾、ａ2＾、・・・、ａ12＾、は、各式において、それぞれ、Ａ、ａの上に＾が付されているものに対応する。係数行列Ａ＾は、（７）式をＡで偏微分した値が０（零）となるときのＡの値として算出され、以下の（１１）式により算出される。従って、係数行列Ａ＾は、以下の（１２）式のようになる。（１１）式および（１２）式において、Ｔは、転置行列であることを示す（このことは、その他の数式でも同じである）。
【0052】
【数6】
【0053】
（１２）式の右辺のＣ1
TＣ1∈Ｒ3n1×3n1は、以下の（１３）式で表される。
【0054】
【数7】
【0055】
そこで、以下の（１４）式により定義される行列Ｚ4を用いて、（１２）式の右辺のＣ1
TＣ1を以下の（１５）式のように表すこととする。（１５）式において、行列０4は、４行４列の０行列（全ての成分が０（零）の行列）である。
【0056】
【数8】
【0057】
また、（１２）式の左辺の係数行列Ａ＾Tを、以下の（１６）式のように３つのブロックに分割した形で表現することとする。
【0058】
【数9】
【0059】
同様に、（１２）式の右辺のＣ1
TＣ2を、以下の（１７）式のように３つのブロックに分割した形で表現することとする。
【0060】
【数10】
【0061】
そうすると、（１５）式の関係から、（１１）式は、以下の（１８）式～（２０）式に示す３つの方程式で表される。
【0062】
【数11】
【0063】
（１８）式～（２０）式を式変形し、Σ（積算記号）を用いない形に書き直すと、以下の（２１）式～（２３）式が得られる。
【0064】
【数12】
【0065】
ここで、部位の平行移動が行われても、当該部位には応力は作用しない。一方、部位に応力が作用すると、当該部位の線形変換（拡縮、回転、およびせん断のうちの少なくとも１つ）が行われ。そこで、本実施形態では、係数行列Ａ＾の成分のうち、応力に関係するアフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を、以下の（２４）式のように行列Ｍ∈Ｒ9×1と定義する。ここで、アフィン係数ａ4＾、ａ8＾、ａ12＾は、アフィン変換行列の成分のうち、平行移動を行うための成分である。アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾は、線形変換を行うための成分である。本実施形態では、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾が、物体（本実施形態では台車枠１６）の応力に対する影響因子である場合を例示する。また、本実施形態では、台車枠１６の部位ｉの座標（位置）が、第２物理量である場合を例示する。また、本実施形態では、アフィン変換前の台車枠１６の部位ｉの座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）が、変換前の第２物理量であり、アフィン変換後の台車枠１６の部位ｉの座標（ｘ2
(i)，ｙ2
(i)，ｚ2
(i)）が、変換後の第２物理量である場合を例示する。
【0066】
【数13】
【0067】
（２１）式～（２３）式より、行列Ｍは、以下の（２５）式で表される。本実施形態では、（３）式や（２５）式が、変換前の第２物理量と、当該変換前の第２物理量に乗算される係数と、に基づいて、変換後の前記第２物理量を算出する計算式である場合を例示する。
【0068】
【数14】
【0069】
前述したように本実施形態では、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾が、線形回帰式における複数の説明変数であり、台車枠１６の最大主応力が線形回帰式における目的変数である場合を例示する。また、説明の都合上、台車枠１６の部位ｊにおける最大主応力をσmax,jと表記する。また、台車枠１６の部位ｊにおける最大主応力をσmax,jを算出するための線形回帰式の回帰係数をｗ0,j、ｗ1,j、・・・、ｗ9,jと表記する。また、台車枠１６の部位ｊに応じて定められたｎ1個の部位ｉ（ｉ＝１，２，・・・ｎ1）の座標（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）、（ｘ2
(i)，ｙ2
(i)，ｚ2
(i)）に基づいて算出されるアフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を、それぞれａ1
(j)＾、ａ2
(j)＾、ａ3
(j)＾、ａ5
(j)＾、ａ6
(j)＾、ａ7
(j)＾、ａ9
(j)＾、ａ10
(j)＾、ａ11
(j)＾と表記する。各変数において上付きの（ｊ）は、部位ｊの値であることを表す記号である。
【0070】
本実施形態では、台車枠１６の部位ｊにおける最大主応力σmax,jが、以下の（２６）式の線形回帰式で表される場合を例示する。なお、各変数において下付きの“ｊ”は、部位ｊの値であることを表す記号である。台車枠１６の部位ｊは、例えば、メッシュの中から選ばれる。また、部位ｊの位置は、台車枠１６の応力（本実施形態では最大主応力）の算出対象の位置に対応する。線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを算出する際には、部位ｊの位置は、台車枠１６の可及的に広範囲の領域に分布するようにするのが好ましい。線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jの計算精度が高まるからである。
【0071】
【数15】
【0072】
なお、回帰係数算出部４２０により、線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを算出した後、後述する応力算出部４４０により、台車枠１６の部位ｊにおける最大主応力σmax,jを推定する際には、部位ｊの数は、最大主応力σmax,jを算出する部位の数に応じた数になるので、１つでもよいし、複数でもよい。また、台車枠１６の部位ｊにおける最大主応力σmax,jを推定する際に用いる部位ｊは、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを算出する際に用いる部位ｉに含まれていてもよいし、含まれていなくてもよい。
【0073】
ここで、説明の都合上（目的変数および説明変数が学習時および推定時の何れにおいて用いられるものであるのかを区別するために）、教師データを識別する記号をｐと表記し、教師データｐに含まれる台車枠１６の部位ｊにおける最大主応力をσmax,p,jと表記する。また、教師データｐに含まれる台車枠１６の部位ｊにおける最大主応力をσmax,p,jを、以下の（２７）式のようにψp,jと表記する。また、教師データｐに含まれる（２６）式の右辺の説明変数の成分ａ1
(j)＾、ａ2
(j)＾、ａ3
(j)＾、ａ5
(j)＾、ａ6
(j)＾、ａ7
(j)＾、ａ9
(j)＾、ａ10
(j)＾、ａ11
(j)＾を、それぞれ、ａ1
(p,j)＾、ａ2
(p,j)＾、ａ3
(p,j)＾、ａ5
(p,j)＾、ａ6
(p,j)＾、ａ7
(p,j)＾、ａ9
(p,j)＾、ａ10
(p,j)＾、ａ11
(p,j)＾と表記する。また、教師データｐに含まれる（２６）式の右辺の説明変数行列［ａ1
(j)＾ ａ2
(j)＾ ａ3
(j)＾ ａ5
(j)＾ ａ6
(j)＾ ａ7
(j)＾ ａ9
(j)＾ ａ10
(j)＾ ａ11
(j)＾］Tを、以下の（２８）式の左辺のようにλp,jと表記する。なお、各変数において下付きの"ｐ"は、教師データｐの値であることを表す記号である。また、各変数において上付きの（ｐ，ｊ）は、教師データｐに含まれる部位ｊの値であることを表す記号である。
【0074】
また、教師データｐに含まれる複数の説明変数であるアフィン係数ａ1
(pj)＾、ａ2
(p,j)＾、ａ3
(p,j)＾、ａ5
(p,j)＾、ａ6
(p,j)＾、ａ7
(p,j)＾、ａ9
(p,j)＾、ａ10
(p,j)＾、ａ11
(p,j)＾を、それぞれ、以下の（２８）式のようにλp,j,1、λp,j,2、λp,j,3、λp,j,4、λp,j,5、λp,j,6、λp,j,7、λp,j,8、λp,j,9、と表記する。また、（２６）式の右辺の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを成分として含む係数行列を、以下の（２９）式のようにＷjと表記する。
【0075】
【数16】
【0076】
係数行列Ｗは、重回帰分析における教師あり学習を行うことにより算出される。ここでは、最小二乗法により係数行列Ｗjを算出する場合を例示する。この場合、係数行列Ｗjは、以下の（３０）式により算出される。なお、（３０）式は、（３０）式の［］内の値が最小になるときの係数行列Ｗjを算出することを示す数式である。
【0077】
【数17】
【0078】
ここで、（３０）式の［］内の値が最小になるときの係数行列Ｗjを係数行列Ｗj＾と表記し、係数行列Ｗj＾の成分である回帰係数をｗ0,j＾、ｗ1,j＾、・・・、ｗ9,j＾と表記する（以下の（３１）式を参照）。なお、Ｗj＾、ｗ0,j＾、ｗ1,j＾、・・・、ｗ9,j＾は、各式において、それぞれ、Ｗ、ｗの上に＾が付されているものに対応する。係数行列Ｗj＾は、（３０）式をＷで偏微分した値が０（零）となるときのＷの値として算出され、以下の（３２）式の関係を満たす。ここで、ｎ2は、教師データの数である。ｎ2の数を多くすれば、線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jの計算精度は高まるが計算負荷も高まる。このような観点からｎ2は予め設定される。
【0079】
【数18】
【0080】
（３２）式を、ｎ2個の数式に分けて記載すると、以下の（３３）式のようになる。なお、表記の都合上、（３３）式では、ｐ＝２～ｎ2－１までの数式を示していない。
【0081】
【数19】
【0082】
（３３）式の行列内の演算を展開すると、以下の（３４）式が得られる。なお、表記の都合上、（３４）式においても、ｐ＝２～ｎ2－１までの数式を示していない。
【0083】
【数20】
【0084】
（３４）式に示すｎ2個の数式の行列和をとると、以下の（３５）式が得られる。さらに、以下の（３５）式を式展開すると、以下の（３６）式が得られる。
【0085】
【数21】
【0086】
ここで、（３６）式の右辺の１つ目の行列、２つ名の行列を、それぞれ、以下の（３７）式、（３８）式のように行列Ψj、Λjと表記する。そうすると、（３６）式は、以下の（３９）式のように表される。したがって、係数行列Ｗj＾は、以下の（４０）式または（４１）式で表される。
【0087】
【数22】
【0088】
本実施形態では、回帰係数算出部４２０が（４０）式または（４１）式により係数行列Ｗj＾を算出する場合を例示する。係数行列Ｗj＾を算出する方法の具体例を以下に説明する。
まず、回帰係数算出部４２０は、台車枠１６の部位ｊとして少なくとも１つの部位ｊを設定する。また、回帰係数算出部４２０は、１つの部位ｊに対してそれぞれｎ1個の部位ｉを設定する。また、回帰係数算出部４２０は、台車枠１６の各メッシュの初期座標（ｘ0，ｙ0，ｚ0）を設定する。
まず、回帰係数算出部４２０は、台車枠１６の部位ｊとして少なくとも１つの部位ｊを設定する。また、回帰係数算出部４２０は、１つの部位ｊに対してそれぞれｎ1個の部位ｉを設定する。また、回帰係数算出部４２０は、台車枠１６の各メッシュの初期座標（ｘ0，ｙ0，ｚ0）を設定する。
【0089】
本実施形態では、回帰係数算出部４２０が、（２）式に基づいて算出される外力ベクトル｛ｆ｝が与えられた場合の（１）式の運動方程式を、有限要素法により解くことで、台車枠１６の解析モデルにおける各メッシュの節点の位置での変位ベクトル[ｕ]を算出することにより、台車枠１６の部位ｉ、ｊの座標を算出する場合を例示する。
【0090】
具体的に回帰係数算出部４２０は、加速度センサ２１ａ、２１ｂで測定された軸箱１７ａ、１７ｂの加速度データを用いて、軸箱１７ａ、１７ｂの加速度の時間に関する１階積分、２階積分を行うことにより、（２）式の左辺の速度ベクトル｛ｕ0・｝、変位ベクトル｛ｕ0｝をそれぞれ算出する。ｕ0・は、（２）式においてｕ0の上に・が付されているものに対応する。回帰係数算出部４２０は、モノリンク１８ａ、１８ｂの粘性減衰係数、軸バネ１９ａ、１９ｂの粘性減衰係数から、モノリンク１８ａ、１８ｂの粘性行列［Ｃbc］、軸バネ１９ａ、１９ｂの粘性行列［Ｃbc］をそれぞれ算出する。また、回帰係数算出部４２０は、モノリンク１８ａ、１８ｂの剛性、軸バネ１９ａ、１９ｂの剛性から、モノリンク１８ａ、１８ｂの剛性行列［Ｋbc］、軸バネ１９ａ、１９ｂの剛性行列［Ｋbc］をそれぞれ算出する。
【0091】
回帰係数算出部４２０は、このようにして得られた情報を（２）式に与えることにより、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力と、軸バネ１９ａ、１９ｂの着力点３２ａ、３２ｂに作用する外力とを算出する。そして、回帰係数算出部４２０は、台車枠１６の着力点３１ａ、３１ｂに作用する外力、着力点３２ａ、３２ｂに作用する外力を、それぞれ、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力の反作用力、軸バネ１９ａ、１９ａの着力点３２ａ、３２ｂに作用する外力の反作用力として算出する。
【0092】
回帰係数算出部４２０は、（１）式の外力ベクトル｛ｆ｝の成分のうち、着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力の成分には、このようにして算出した値を与え、その他の成分には、０（零）を与えることにより、（１）式の外力ベクトル｛ｆ｝を算出する。そして、回帰係数算出部４２０は、外力ベクトル｛ｆ｝が与えられた（１）式の運動方程式を解くことにより、台車枠１６の解析モデルにおける各メッシュの節点の位置での変位ベクトル｛ｕ]を算出する。なお、台車枠１６の質量行列［Ｍ］、粘性行列［Ｃ］、および剛性行列［Ｋ］は、有限要素法による運動方程式の定式化に従って算出される。また、台車枠１６の解析モデルにおける各メッシュの節点の位置での変位ベクトル｛ｕ]を算出する手法自体は公知の技術で実現され、前述した手法に限定されない。
【0093】
また、本実施形態では、回帰係数算出部４２０は、加速度センサ２２ａ、２２ｂで測定された台車枠１６の加速度データを用いて、台車枠１６の加速度の時間に関する２階積分を行うことにより、加速度センサ２２ａ、２２ｂの取り付け位置における変位ベクトル（ｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分）を算出する。そして、回帰係数算出部４２０は、台車枠１６の加速度データを用いて算出された、加速度センサ２２ａ、２２ｂの取り付け位置での変位ベクトルと、（１）式の運動方程式を解くことにより算出された、当該取り付け位置での変位ベクトルと、の差に基づいて、（１）式の運動方程式を解くことにより算出された各メッシュの節点の位置における変位ベクトルを補正する。例えば、回帰係数算出部４２０は、（１）式の運動方程式を解くことにより算出された変位ベクトルが、台車枠１６の加速度データを用いて算出された変位ベクトルに一致するように、各メッシュの節点の位置での変位ベクトルに対する補正量を算出する。そして、回帰係数算出部４２０は、（１）式の運動方程式を解くことにより算出された変位ベクトルを、当該補正量を用いて補正する。なお、変位ベクトルの補正は行われなくてもよい。
【0094】
回帰係数算出部４２０は、鉄道車両が走行しているときに、台車枠１６の各メッシュの節点の位置における変位ベクトル｛ｕ]を、所定の時間隔Δｔ１が経過するたびに算出する。そして、回帰係数算出部４２０は、台車枠１６の各メッシュの初期座標（ｘ0，ｙ0，ｚ0）および変位ベクトル｛ｕ]に基づいて、台車枠１６の各メッシュの節点の位置における座標（ｘ，ｙ，ｚ）を、所定の時間隔Δｔ１の各学習時刻ｔ１において算出する。このようにして、台車枠１６の各メッシュの節点の位置における座標（ｘ，ｙ，ｚ）が各学習時刻ｔ１において算出される。
【0095】
また、回帰係数算出部４２０は、台車枠１６の各メッシュの節点の位置における応力（本実施形態では最大主応力）を各学習時刻ｔ１において算出する。台車枠１６の各メッシュの節点の位置における応力の算出は、特許文献１、２に記載の技術等、公知の技術を用いて行えばよいので、ここでは、その詳細な説明を省略する。
【0096】
回帰係数算出部４２０は、台車枠１６の各メッシュの節点の位置における座標のうち、部位ｊに応じて定められたｎ1個の部位ｉの座標を抽出することと、当該部位ｊに対応する座標における最大主応力σmax,jを抽出することと、を学習時刻ｔ１において行う。回帰係数算出部４２０は、学習時刻ｔ１－Δｔ１におけるｎ1個の部位ｉの座標を（ｘ1
(i)，ｙ1
(i)，ｚ1
(i)）とし、学習時刻ｔ１におけるｎ1個の部位ｉの座標を（ｘ2
(i)，ｙ2
(i)，ｚ2
(i)）として、（２５）式の計算を行うことにより、（２４）式に示す行列Ｍを、学習時刻ｔ１における部位ｊの行列Ｍとして算出する。以上のようにして、学習時刻ｔ１において、部位ｊの最大主応力σmax,p,jと、部位ｊの行列Ｍと、が得られる。回帰係数算出部４２０は、このようにして得られる部位ｊの最大主応力σmax,p,jおよび行列Ｍを、学習時刻ｔ１における部位ｊの教師データｐとする。なお、（２６）式に示すアフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾は、それぞれ、（２８）式に示すａ1
(p,j)＾、ａ2
(p,j)＾、ａ3
(p,j)＾、ａ5
(p,j)＾、ａ6
(p,j)＾、ａ7
(p,j)＾、ａ9
(p,j)＾、ａ10
(p,j)＾、ａ11
(p,j)＾として算出される。
【0097】
回帰係数算出部４２０は、以上の教師データｐの作成を、１つの学習時刻ｔ１において、部位ｊのそれぞれに対して行う。したがって、１つの学習時刻ｔ１において、部位ｊの数の教師データが算出される。また、回帰係数算出部４２０は、各学習時刻ｔ１のそれぞれにおいて、以上の教師データｐを算出することで、１つの部位ｊについて合計でｎ2個の教師データｐを作成する。
【0098】
そして、回帰係数算出部４２０は、ｎ2個の教師データｐを用いて（４０）式または（４１）式の計算を行うことにより係数行列Ｗj＾を算出する。回帰係数記憶部４３０は、以上のようにして回帰係数算出部４２０により算出された係数行列Ｗj＾を記憶する。このようにして（２６）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jが得られる。また、台車枠１６の応力の算出対象の位置に対応する部位ｊが複数ある場合、回帰係数算出部４２０は、複数の部位ｊのそれぞれについて、前述したようにして（２６）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを算出して、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jを成分として含む係数行列Ｗj＾を回帰係数記憶部４３０に記憶させる。
【0099】
なお、行列Ｍおよび係数行列Ｗj＾を算出する手法は、前述した手法に限定されず、重回帰分析における教師あり学習を行うための公知の種々の手法でもよい。
【0100】
<<応力算出部４４０、ステップＳ６０２>>
応力算出部４４０は、回帰係数記憶部４３０により係数行列Ｗj＾が記憶された後に起動する。応力算出部４４０は、データ取得部４１０により取得されたデータと、線形回帰式と、に基づいて、台車枠１６の応力を算出する。本実施形態では、応力算出部４４０が、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jが与えられた（２６）式の線形回帰式に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを算出（推定）する場合を例示する。
応力算出部４４０は、回帰係数記憶部４３０により係数行列Ｗj＾が記憶された後に起動する。応力算出部４４０は、データ取得部４１０により取得されたデータと、線形回帰式と、に基づいて、台車枠１６の応力を算出する。本実施形態では、応力算出部４４０が、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jが与えられた（２６）式の線形回帰式に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを算出（推定）する場合を例示する。
【0101】
応力算出部４４０は、最大主応力σmax,jの算出対象となる部位ｊとして１つ以上の部位ｊを設定する。また、応力算出部４４０は、１つの部位ｊに対してそれぞれｎ1個の部位ｉを設定する。また、応力算出部４４０は、台車枠１６の各メッシュの初期座標（ｘ0，ｙ0，ｚ0）を設定する。
【0102】
そして、応力算出部４４０は、<<回帰係数算出部４２０、回帰係数記憶部４３０、ステップＳ５０２、Ｓ５０３>>の欄で説明した所定の時間隔Δｔ１の学習時刻ｔ１を、所定の時間隔Δｔ２の推定時刻ｔ２に置き替えた処理を行って、推定時刻ｔ２における部位ｊの行列Ｍを算出する（（２４）式を参照）。最大主応力σmax,jの算出対象となる部位ｊが複数ある場合、応力算出部４４０は、各推定時刻ｔ２において、各部位ｊの行列Ｍをそれぞれ算出する。
【0103】
そして、応力算出部４４０は、各推定時刻ｔ２における部位ｊの行列Ｍと、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ9,jが与えられた（２６）式の線形回帰式と、に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを、各推定時刻ｔ２において算出する。最大主応力σmax,jの算出対象となる部位ｊが複数ある場合、応力算出部４４０は、各推定時刻ｔ２において、各部位ｊにおける最大主応力σmax,jを算出する。
【0104】
<<出力部４５０、ステップＳ６０３>>
出力部４５０は、応力算出部４４０により算出された台車枠１６の部位ｊにおける最大主応力σmax,jを示す情報を出力する。出力部４５０による情報の出力形態は限定されない。例えば、出力部４５０による情報の出力形態として、コンピュータディスプレイへの表示と、外部装置への送信と、処理装置４００の内部または外部の記憶媒体への記憶と、のうちの少なくとも１つが採用される。
出力部４５０は、応力算出部４４０により算出された台車枠１６の部位ｊにおける最大主応力σmax,jを示す情報を出力する。出力部４５０による情報の出力形態は限定されない。例えば、出力部４５０による情報の出力形態として、コンピュータディスプレイへの表示と、外部装置への送信と、処理装置４００の内部または外部の記憶媒体への記憶と、のうちの少なくとも１つが採用される。
【0105】
＜計算例＞
次に、本実施形態の計算例を説明する。図７は、特許文献１に記載の手法で算出された台車枠１６の部位ｊの最大主応力と時間との関係の一例を示す図である。図７では、４／３６００秒の時間隔で各時刻における最大主応力を算出した結果を示す。図７における前半の６００個の各時刻において、本実施形態で説明したようにして部位ｊの行列Ｍ（アフィン係数ａ1 (p,j)＾、ａ2 (p,j)＾、ａ3 (p,j)＾、ａ5 (p,j)＾、ａ6 (p,j)＾、ａ7 (p,j)＾、ａ9 (p,j)＾、ａ10 (p,j)＾、ａ11 (p,j)＾）を算出した。これら６００個の各時刻における行列Ｍおよび最大主応力を教師データとして、本実施形態で説明したようにして係数行列Ｗj＾を算出した。そして、図７における後半の３００１個の各時刻において、特許文献１に記載の手法と本実施形態の手法とのそれぞれの手法で、台車枠１６の部位ｊの最大主応力を算出した。また、前半の６００個の各時刻においても、前述したようにして算出した係数行列Ｗj＾を用いて本実施形態の手法で台車枠１６の部位ｊの最大主応力を算出した。なお、最大主応力の算出対象の部位ｊは、いずれも同じ部位である。
次に、本実施形態の計算例を説明する。図７は、特許文献１に記載の手法で算出された台車枠１６の部位ｊの最大主応力と時間との関係の一例を示す図である。図７では、４／３６００秒の時間隔で各時刻における最大主応力を算出した結果を示す。図７における前半の６００個の各時刻において、本実施形態で説明したようにして部位ｊの行列Ｍ（アフィン係数ａ1 (p,j)＾、ａ2 (p,j)＾、ａ3 (p,j)＾、ａ5 (p,j)＾、ａ6 (p,j)＾、ａ7 (p,j)＾、ａ9 (p,j)＾、ａ10 (p,j)＾、ａ11 (p,j)＾）を算出した。これら６００個の各時刻における行列Ｍおよび最大主応力を教師データとして、本実施形態で説明したようにして係数行列Ｗj＾を算出した。そして、図７における後半の３００１個の各時刻において、特許文献１に記載の手法と本実施形態の手法とのそれぞれの手法で、台車枠１６の部位ｊの最大主応力を算出した。また、前半の６００個の各時刻においても、前述したようにして算出した係数行列Ｗj＾を用いて本実施形態の手法で台車枠１６の部位ｊの最大主応力を算出した。なお、最大主応力の算出対象の部位ｊは、いずれも同じ部位である。
【0106】
図８にその結果を示す。図８（ａ）および図８（ｂ）において、真値は、特許文献１に記載の手法で算出された台車枠１６の部位ｊの最大主応力である。推定値は、本実施形態の手法で算出された台車枠１６の当該部位ｊの最大主応力である。同時刻における真値および推定値を１つの点としてプロットすることにより、図８（ａ）および図８（ｂ）が得られる。
【0107】
図８（ａ）は、前述した前半の６００個の各時刻における、台車枠１６の最大主応力の真値と推定値との関係の一例を示す図である。図８（ｂ）は、前述した後半の３００１個の各時刻における、台車枠１６の最大主応力の真値と推定値との関係の一例を示す図である。図８（ａ）においては、決定係数Ｒ2は、０．９９８０３３７２９４となり、図８（ｂ）においては、決定係数Ｒ2は、０．９９８０９５８０５９となり、本実施形態の手法では、特許文献１に記載の手法と同程度の精度で、台車枠１６の最大主応力を算出することができることが分かる。また、ＣＰＵの処理時間は、特許文献１に記載の手法よりも本実施形態の手法の方が、２２倍以上速く、また、台車枠１６の最大主応力の予測に要する時間は、特許文献１に記載の手法よりも本実施形態の手法の方が、１０００倍以上速くなった。したがって、本実施形態の手法では、特許文献１に記載の手法に比べ、計算時間を大幅に短縮することができることが分かる。
【0108】
＜まとめ＞
以上のように本実施形態では、処理装置４００は、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データを取得する。処理装置４００は、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データに基づいて、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を算出する。処理装置４００は、台車枠１６の最大主応力σmax,jを目的変数として含み、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を複数の説明変数として含む線形回帰式に、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を与えることにより、台車枠１６の最大主応力σmax,jを算出する。したがって、特許文献１に記載のように、アフィン行列の特異値分解を行ったり、特許文献２に記載のように、モード座標系から物理座標系への再変換を行ったり、特許文献１、２に記載のように、変位から歪みを算出し歪みから応力を算出したりしなくても、複数の説明変数から物体（本実施形態では台車枠１６）の応力を直接的に算出することができる。よって、物体の応力を推定するに際し、計算速度を向上させることができる。
以上のように本実施形態では、処理装置４００は、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データを取得する。処理装置４００は、台車枠１６の加速度データおよび軸箱１７ａ、１７ｂの加速度データに基づいて、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を算出する。処理装置４００は、台車枠１６の最大主応力σmax,jを目的変数として含み、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を複数の説明変数として含む線形回帰式に、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾を与えることにより、台車枠１６の最大主応力σmax,jを算出する。したがって、特許文献１に記載のように、アフィン行列の特異値分解を行ったり、特許文献２に記載のように、モード座標系から物理座標系への再変換を行ったり、特許文献１、２に記載のように、変位から歪みを算出し歪みから応力を算出したりしなくても、複数の説明変数から物体（本実施形態では台車枠１６）の応力を直接的に算出することができる。よって、物体の応力を推定するに際し、計算速度を向上させることができる。
【0109】
＜変形例＞
本実施形態では、台車枠１６の振動（運動）を表す運動方程式を用いて、台車枠１６の変位ベクトル[ｕ]を算出し、台車枠１６の変位ベクトル[ｕ]に基づいて、台車枠１６の各部位ｊの座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）、（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）を算出する場合を例示した。しかしながら、必ずしもこのようにする必要はない。例えば、台車枠１６の各部位ｉの変位または変位を算出可能な物理量（例えば、加速度および速度）を測定するセンサ（例えば、歪みゲージおよび加速度センサ）の測定値に基づいて、当該部位ｉの変位を得て、当該部位ｉの座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）、（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）を算出してもよい。
本実施形態では、台車枠１６の振動（運動）を表す運動方程式を用いて、台車枠１６の変位ベクトル[ｕ]を算出し、台車枠１６の変位ベクトル[ｕ]に基づいて、台車枠１６の各部位ｊの座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）、（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）を算出する場合を例示した。しかしながら、必ずしもこのようにする必要はない。例えば、台車枠１６の各部位ｉの変位または変位を算出可能な物理量（例えば、加速度および速度）を測定するセンサ（例えば、歪みゲージおよび加速度センサ）の測定値に基づいて、当該部位ｉの変位を得て、当該部位ｉの座標（ｘ1 (i)，ｙ1 (i)，ｚ1 (i)）、（ｘ2 (i)，ｙ2 (i)，ｚ2 (i)）を算出してもよい。
【0110】
また、本実施形態では、アフィン変換行列を用いる場合を例示した。しかしながら、必ずしもアフィン変換行列を用いる必要はない。例えば、線形変換行列（アフィン変換行列から平行移動を行うための平行移動を行うための成分を除いた行列）を用いてもよい。
【0111】
また、本実施形態では、目的変数が最大主応力である場合を例示した。しかしながら、目的変数となる応力は、最大主応力に限定されない。目的変数は、例えば、中間主応力でも、最小主応力でも、せん断応力であってもよい。
【0112】
また、本実施形態では、学習を行う機能（図５のフローチャートによる処理）と、推定を行う機能（図６のフローチャートによる処理）とが、１つの装置（処理装置４００）で実現される場合を例示した。しかしながら、必ずしもこのようにする必要はない。例えば、推定を行う機能を有する処理装置（推定装置）が、学習を行う機能を有する装置とは別の装置であってもよい。
【0113】
（第２実施形態）
次に、第２実施形態を説明する。第１実施形態では、複数の説明変数に、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾が含まれる場合を例示した。これに対し、本実施形態では、複数の説明変数に、モード座標系における変位および速度のうちの少なくとも一方が含まれる場合について説明する。このように本実施形態と第１実施形態とは、説明変数が異なることによる構成および処理が主として異なる。したがって、本実施形態の説明において、第１実施形態と同一の部分については、図１～図８に付した符号と同一の符号を付す等して詳細な説明を省略する。
次に、第２実施形態を説明する。第１実施形態では、複数の説明変数に、アフィン係数ａ1＾、ａ2＾、ａ3＾、ａ5＾、ａ6＾、ａ7＾、ａ9＾、ａ10＾、ａ11＾が含まれる場合を例示した。これに対し、本実施形態では、複数の説明変数に、モード座標系における変位および速度のうちの少なくとも一方が含まれる場合について説明する。このように本実施形態と第１実施形態とは、説明変数が異なることによる構成および処理が主として異なる。したがって、本実施形態の説明において、第１実施形態と同一の部分については、図１～図８に付した符号と同一の符号を付す等して詳細な説明を省略する。
【0114】
＜処理装置４００、処理方法＞
処理装置４００のハードウェアは、例えば、第１実施形態で説明したハードウェアで実現される。また、処理装置４００の機能的な構成は、例えば、第１実施形態で説明した図４に示す各ブロック（データ取得部４１０、回帰係数算出部４２０、回帰係数記憶部４３０、応力算出部４４０、および出力部４５０）が有する機能の少なくとも一部が、第１実施形態で説明した機能に対し変更されたもので実現される。また、処理方法は、例えば、図５および図６に示す各ステップで行われる処理の少なくとも一部が、第１実施形態で説明した処理に対し変更されたもので実現される。そこで、図４～図６に示す符号を用いて、本実施形態の処理装置４００および処理方法の一例を説明する。
処理装置４００のハードウェアは、例えば、第１実施形態で説明したハードウェアで実現される。また、処理装置４００の機能的な構成は、例えば、第１実施形態で説明した図４に示す各ブロック（データ取得部４１０、回帰係数算出部４２０、回帰係数記憶部４３０、応力算出部４４０、および出力部４５０）が有する機能の少なくとも一部が、第１実施形態で説明した機能に対し変更されたもので実現される。また、処理方法は、例えば、図５および図６に示す各ステップで行われる処理の少なくとも一部が、第１実施形態で説明した処理に対し変更されたもので実現される。そこで、図４～図６に示す符号を用いて、本実施形態の処理装置４００および処理方法の一例を説明する。
【0115】
＜モード座標系における変位および速度の算出＞
本実施形態では、以下のようにしてモード座標系における物体の変位ベクトル｛ξ｝および速度ベクトル｛ξ・｝を算出する場合を例示する。なお、ここでは、特許文献２に記載されている事項についての詳細な説明を必要に応じて省略する。また、モード座標系における物体の変位ベクトル｛ξ｝（∈Ｒn）および速度ベクトル｛ξ・｝（∈Ｒn）の算出は、例えば、特許文献２に記載の手法を用いてもよい。
本実施形態では、以下のようにしてモード座標系における物体の変位ベクトル｛ξ｝および速度ベクトル｛ξ・｝を算出する場合を例示する。なお、ここでは、特許文献２に記載されている事項についての詳細な説明を必要に応じて省略する。また、モード座標系における物体の変位ベクトル｛ξ｝（∈Ｒn）および速度ベクトル｛ξ・｝（∈Ｒn）の算出は、例えば、特許文献２に記載の手法を用いてもよい。
【0116】
（１）式の運動方程式に固有値解析を適用することにより、固有振動数ωおよび固有ベクトル｛φ｝が算出され、固有ベクトル｛φ｝＝｛φ(1)｝、・・・、｛φ(n)｝で構成されるモード行列［φ］（∈Ｒ3l×n）が算出される。尚、固有ベクトルは、固有モードベクトルとも称される。ここで、ｎ（∈Ｎ）は、モード数である。固有ベクトル｛φ(1)｝、・・・、｛φ(n)｝を基底とする変位ベクトルの座標系のことをモード座標系と称する。
【0117】
モード座標系における台車枠１６の質量行列［Ｍξ］（∈Ｒn×n）、モード座標系における台車枠１６の粘性行列［Ｃξ］（∈Ｒn×n）、モード座標系における台車枠１６の剛性行列［Ｋξ］（∈Ｒn×n）は、それぞれ、以下の（４２）式、（４３）式、（４４）式で表される。
【0118】
【数23】
【0119】
また、モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］は、それぞれ、以下の（４５）式、（４６）式、（４７）式のように、対角行列である。
【0120】
【数24】
【0121】
ここで、（１）、・・・、（ｎ）は、それぞれ、一次の固有振動モード、・・・、ｎ次の固有振動モードに対応する成分であることを示す（このことは、以降の式でも同じである）。
モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］は、対角行列である。従って、各固有振動モードは、相互に独立したものとして扱うことができる。よって、計算時間を短くすることができる。
モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］は、対角行列である。従って、各固有振動モードは、相互に独立したものとして扱うことができる。よって、計算時間を短くすることができる。
【0122】
物理座標系における台車枠１６の変位ベクトル｛ｕ｝は、以下の（４８）式のように、モード座標系における変位ベクトル｛ξ｝（∈Ｒn）に変換される。また、物理座標系における台車枠１６の速度ベクトル｛ｕ・｝は、以下の（４９）式のように、モード座標系における速度ベクトル｛ξ・｝（∈Ｒn）に変換される。また、物理座標系における台車枠１６の外力ベクトル｛ｆ｝は、以下の（５０）式のように、モード座標系における外力ベクトル｛ｆξ｝（∈Ｒn）に変換される。
【0123】
【数25】
【0124】
尚、ｕ・は、（１）式および（４９）式においてｕの上に・が付されているものに対応する。また、ξ・は、（４９）式においてξの上に・が付されているものに対応する（このことは、以降の式でも同じである）。
【0125】
本実施形態では、物理座標系における変位および速度が、変換前の第２物理量である場合を例示する。また、本実施形態では、モード座標系における変位および速度が、変換後の第２物理量、および物体（本実施形態では台車枠１６）の応力に対する影響因子である場合を例示する。また、本実施形態では、（４８）式および（４９）式が、変換前の第２物理量と、当該変換前の第２物理量に乗算される係数と、に基づいて、変換後の前記第２物理量を算出する計算式である場合を例示する。
【0126】
台車枠１６の振動を表す物理座標系における運動方程式（（１）式）は、以下の（５１）式のように、モード座標系における運動方程式で表現される。
【0127】
【数26】
【0128】
また、状態変数Ξ（∈Ｒ2n）を以下の（５２）式のように定義する。また、（５１）式を以下の（５３）式の状態方程式で記述する。尚、特許文献１に記載されているように、（５３）式の右辺にシステムノイズが加算されてもよいが、ここでは、説明を簡単にするために、システムノイズの表記を省略する。
【0129】
【数27】
【0130】
ここで、状態遷移行列Ｓ（∈Ｒ2n×2n）は、以下の（５４）式で表される。
【0131】
【数28】
【0132】
また、ベクトルＦ（∈Ｒn）は、モード座標系における台車枠１６の外力ベクトル｛ｆξ｝を格納するベクトルであり、以下の（５５）式で表される。
【0133】
【数29】
【0134】
また、行列Ｇ（∈Ｒ2n×n）は、以下の（５６）式で表される。
【0135】
【数30】
【0136】
本実施形態では、（５３）式を離散化した式として以下の（５７）式を用いる場合を例示する。ここで、ｋは、離散化された時刻を識別する変数である。また、（５７）式の右辺の係数Ｂ、Γは、それぞれ以下の（５８）式、（５９）式で表される。
【0137】
【数31】
【0138】
ここで、ｅは、ネイピア数である。Ｔsは、サンプリング時間である。サンプリング時間Ｔsは、例えば、時刻ｋ、ｋ－１の時間隔をｑ個に等分した場合の各部分の時間である。また、初期時刻（ｋ＝０）における状態変数Ξ0には、所定値（例えば、０（零））が与えられる。したがって、例えば、ｋが所定値未満の状態変数Ξkを、信頼性が低いものとして採用せずに、時刻ｋが当該所定値以上の状態変数Ξkを採用してもよい。
【0139】
（５７）式を解くことにより各時刻ｋにおける状態変数Ξkが算出され、モード座標系における変位ベクトル[ξ]および速度ベクトル[ξ・]が算出される。なお、モード座標系における変位ベクトル[ξ]および速度ベクトル[ξ・]は、各メッシュのそれぞれにおいて算出される。
【0140】
<<データ取得部４１０、ステップＳ５０１、Ｓ６０１>>
データ取得部４１０は、処理装置４００における計算のために予め処理装置４００が取得しておく必要がある各種の情報を取得する。本実施形態においても、ステップＳ５０１、６０１において、データ取得部４１０は、第１実施形態で説明したデータを取得する。本実施形態では、軸箱１７ａ、１７ｂの加速度データが、物体の振動に応じて変化する第１物理量の測定データである場合を例示する。
データ取得部４１０は、処理装置４００における計算のために予め処理装置４００が取得しておく必要がある各種の情報を取得する。本実施形態においても、ステップＳ５０１、６０１において、データ取得部４１０は、第１実施形態で説明したデータを取得する。本実施形態では、軸箱１７ａ、１７ｂの加速度データが、物体の振動に応じて変化する第１物理量の測定データである場合を例示する。
【0141】
<<回帰係数算出部４２０、回帰係数記憶部４３０、ステップＳ５０２、Ｓ５０３>>
回帰係数算出部４２０は、データ取得部４１０により取得されたデータを用いて、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式の回帰係数を算出する。本実施形態でも第１実施形態と同様に、台車枠１６の応力（目的変数）が、台車枠１６の最大主応力である場合を例示する。また、本実施形態では、複数の説明変数が、状態変数Ξ（＝［ξ・(1) ξ(1) ・・・ξ・(n) ξ(n)］Tである場合を例示する。
回帰係数算出部４２０は、データ取得部４１０により取得されたデータを用いて、台車枠１６の応力を目的変数として含み、台車枠１６の応力に対する影響因子を複数の説明変数として含む線形回帰式の回帰係数を算出する。本実施形態でも第１実施形態と同様に、台車枠１６の応力（目的変数）が、台車枠１６の最大主応力である場合を例示する。また、本実施形態では、複数の説明変数が、状態変数Ξ（＝［ξ・(1) ξ(1) ・・・ξ・(n) ξ(n)］Tである場合を例示する。
【0142】
ここで、台車枠１６の応力の算出対象の位置に対応する部位ｊの状態変数を、Ξj（＝［ξ・j
(1) ξj
(1)・・・ξ・j
(n) ξj
(n)］Tと表記する。本実施形態では、台車枠１６の部位ｊにおける最大主応力σmax,jが、以下の（６０）式の線形回帰式で表される場合を例示する。
【0143】
【数32】
【0144】
説明の都合上、教師データｐに含まれる（６０）式の右辺の説明変数の成分ξ・j
(1)、ξj
(1)、・・・、ξ・j
(n)、ξj
(n)を、それぞれ、ξ・p,j
(1)、ξp,j
(1)、・・・、ξ・p,j
(n)、ξp,j
(n)と表記する。また、第１実施形態との重複する説明を省略するため、ここでは第１実施形態で説明した記号を用いる。
【0145】
そうすると、係数行列Ｗj＾は、以下の（６１）式で表される。また、係数行列Ｗj＾を算出するために用いる教師データｐに含まれる目的変数（台車枠１６の部位ｊにおける最大主応力σmax,p,j）を成分として含む行列Ψjは、第１実施形態で説明した（３７）式で表される。また、係数行列Ｗj＾を算出するために用いる教師データｐに含まれる複数の説明変数（状態変数（ξ・p,j
(1)、ξp,j
(1)、・・・、ξ・p,j
(n)、ξp,j
(n)））を成分として含む行列Λjは、以下の（６２）式で表される（ｐ＝１，２，・・・，ｎ2）。（６１）式、（６２）式は、それぞれ、第１実施形態で説明した（３１）式、（３８）式に対応する。
【0146】
【数33】
【0147】
そうすると、（３７）式、（６２）式で表される行列Ψj、Λjを用いて、第１実施形態で説明したようにして（４０）式または（４１）式の計算を行うことにより、（６１）式に示す係数行列Ｗj＾が算出される。
【0148】
また、複数の説明変数は、モード座標系における変位ベクトル[ξ]および速度ベクトル[ξ・]のいずれか一方のみを含んでいてもよい。例えば、複数の説明変数が、モード座標系における変位ベクトル[ξ]を含み、速度ベクトル[ξ・]を含まない場合、（６０）式に示す台車枠１６の部位ｊにおける最大主応力σmax,jを算出するための線形回帰式は、以下の（６３）式のようになる。また、（３１）式で表される係数行列Ｗj＾は、以下の（６４）式のようになる。また、（６２）式で表される行列Λjは、以下の（６５）式のようになる。また、行列Ψjは、（３７）式になる。
【0149】
【数34】
【0150】
この場合、（３７）式、（６５）式で表される行列Ψj、Λjを用いて、（４０）式または（４１）式により、（６４）式に示す係数行列Ｗj＾が算出される。
本実施形態では、回帰係数算出部４２０が以上のようにして（４０）式または（４１）式の計算を行うことにより係数行列Ｗj＾を算出して回帰係数記憶部４３０に記憶させる場合を例示する。
本実施形態では、回帰係数算出部４２０が以上のようにして（４０）式または（４１）式の計算を行うことにより係数行列Ｗj＾を算出して回帰係数記憶部４３０に記憶させる場合を例示する。
【0151】
係数行列Ｗj＾を算出する方法の具体例を以下に説明する。
まず、回帰係数算出部４２０は、台車枠１６の部位ｊとして少なくとも１つの部位ｊを設定する。
回帰係数算出部４２０は、（１）式の運動方程式に対する固有値解析を行うことにより、固有振動数ωおよび固有ベクトル｛φ｝を導出し、固有ベクトル｛φ｝で構成されるモード行列［φ］を算出する。
そして、回帰係数算出部４２０は、モード行列［φ］と、台車枠１６の質量行列［Ｍ］と、を用いて、（４２）式により、モード座標系における台車枠１６の質量行列［Ｍξ］を導出する。以下の説明では、モード座標系における台車枠１６の質量行列を、必要に応じて、モード質量行列と称する。
まず、回帰係数算出部４２０は、台車枠１６の部位ｊとして少なくとも１つの部位ｊを設定する。
回帰係数算出部４２０は、（１）式の運動方程式に対する固有値解析を行うことにより、固有振動数ωおよび固有ベクトル｛φ｝を導出し、固有ベクトル｛φ｝で構成されるモード行列［φ］を算出する。
そして、回帰係数算出部４２０は、モード行列［φ］と、台車枠１６の質量行列［Ｍ］と、を用いて、（４２）式により、モード座標系における台車枠１６の質量行列［Ｍξ］を導出する。以下の説明では、モード座標系における台車枠１６の質量行列を、必要に応じて、モード質量行列と称する。
【0152】
また、回帰係数算出部４２０は、モード行列［φ］と、台車枠１６の粘性行列［Ｃ］と、を用いて、（４３）式により、モード座標系における台車枠１６の粘性行列［Ｃξ］を算出する。以下の説明では、モード座標系における台車枠１６の粘性行列を、必要に応じて、モード粘性行列と称する。
【0153】
また、回帰係数算出部４２０は、モード行列［φ］と、台車枠１６の剛性行列［Ｋ］と、を用いて、（４４）式により、モード座標系における台車枠１６の剛性行列［Ｋξ］を導出する。以下の説明では、モード座標系における台車枠１６の剛性行列を、必要に応じて、モード剛性行列と称する。
【0154】
また、回帰係数算出部４２０は、第１実施形態で説明したようにして加速度センサ２１ａ、２１ｂで測定された加速度データに基づいて、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ0｝および速度ベクトル｛ｕ0・｝をそれぞれ算出する。回帰係数算出部４２０は、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂの粘性減衰係数から、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂの粘性行列［Ｃbc］をそれぞれ導出する。また、外力導出部４０６は、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂの剛性から、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂの剛性行列［Ｋbc］をそれぞれ導出する。
【0155】
また、回帰係数算出部４２０は、加速度センサ２１ａ、２１ｂで測定された軸箱１７ａ、１７ｂの加速度データを用いて、第１実施形態で説明したようにして（１）式の外力ベクトル｛ｆ｝を算出する。そして、回帰係数算出部４２０は、台車枠１６の外力ベクトル｛ｆ｝と、モード行列［φ］と、を用いて、モード座標系における台車枠１６の外力ベクトル｛ｆξ｝（＝［φ］T｛ｆ｝）を算出する。以下の説明では、モード座標系における台車枠１６の外力ベクトルを、必要に応じて、モード外力ベクトルと称する。
【0156】
回帰係数算出部４２０は、以上のようにして、時刻ｋにおいて、モード質量行列［Ｍξ］、モード粘性行列［Ｃξ］、モード剛性行列［Ｋξ］、モード外力ベクトル｛ｆξ｝を、台車枠１６の解析モデルにおける各メッシュのそれぞれにおいて算出し、これらを用いて、（５７）式により、学習時刻ｔ１（時刻ｋ）における状態変数Ξkを、台車枠１６の解析モデルにおける各メッシュのそれぞれにおいて算出する。そして、回帰係数算出部４２０は、これら台車枠１６の解析モデルにおける各メッシュの状態変数Ξkのうち、部位ｊに対応する状態変数状態変数Ξkを、学習時刻ｔ１における部位ｊの状態変数Ξk（ξ・p,j
(1)、ξp,j
(1)、・・・、ξ・p,j
(n)、ξp,j
(n)）とする。
【0157】
また、回帰係数算出部４２０は、台車枠１６の各メッシュの節点の位置における応力（本実施形態では最大主応力）を各学習時刻ｔ１において算出する。第１実施形態で説明したように、台車枠１６の各メッシュの節点の位置における応力の算出は、特許文献１、２に記載の技術等、公知の技術を用いて行えばよいので、ここでは、その詳細な説明を省略する。そして、回帰係数算出部４２０は、学習時刻ｔ１において算出した、台車枠１６の各メッシュの節点の位置における最大主応力のうち、部位ｊに対応する座標における最大主応力を抽出し、学習時刻ｔ１における部位ｊの最大主応力σmax,p,jとする。
【0158】
回帰係数算出部４２０は、以上のようにして学習時刻ｔ１（時刻ｋ）における部位ｊの状態変数Ξk（ξ・p,j
(1)、ξp,j
(1)、・・・、ξ・p,j
(n)、ξp,j
(n)）および最大主応力σmax,p,jを、学習時刻ｔ１における部位ｊの教師データｐとする。
【0159】
回帰係数算出部４２０は、以上の教師データｐの作成を、１つの学習時刻ｔ１において、部位ｊのそれぞれに対して行う。したがって、１つの学習時刻ｔ１において、部位ｊの数の教師データが算出される。また、回帰係数算出部４２０は、各学習時刻ｔ１のそれぞれにおいて、以上の教師データｐを算出することで、１つの部位ｊについて合計でｎ2個の教師データｐを作成する。
【0160】
そして、回帰係数算出部４２０は、ｎ2個の教師データｐを用いて（４０）式または（４１）式の計算を行うことにより係数行列Ｗj＾を算出する。回帰係数記憶部４３０は、以上のようにして回帰係数算出部４２０により算出された係数行列Ｗj＾を記憶する。このようにして（６０）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jまたは（６３）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jが得られる。また、台車枠１６の応力の算出対象の位置に対応する部位ｊが複数ある場合、回帰係数算出部４２０は、複数の部位ｊのそれぞれについて、前述したようにして（６０）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jまたは（６３）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jを算出して、当該回帰係数を成分として含む係数行列Ｗj＾を回帰係数記憶部４３０に記憶させる。
【0161】
なお、係数行列Ｗj＾を算出する手法は、前述した手法に限定されず、重回帰分析における教師あり学習を行うための公知の種々の手法でもよい。
【0162】
<<応力算出部４４０、ステップＳ６０２>>
応力算出部４４０は、回帰係数記憶部４３０により係数行列Ｗj＾が記憶された後に起動する。応力算出部４４０は、データ取得部４１０により取得されたデータと、線形回帰式と、に基づいて、台車枠１６の応力を算出する。本実施形態では、応力算出部４４０が、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jが与えられた（６０）式の線形回帰式、または、回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jが与えられた（６３）式の線形回帰式に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを算出する場合を例示する。
応力算出部４４０は、回帰係数記憶部４３０により係数行列Ｗj＾が記憶された後に起動する。応力算出部４４０は、データ取得部４１０により取得されたデータと、線形回帰式と、に基づいて、台車枠１６の応力を算出する。本実施形態では、応力算出部４４０が、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jが与えられた（６０）式の線形回帰式、または、回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jが与えられた（６３）式の線形回帰式に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを算出する場合を例示する。
【0163】
応力算出部４４０は、最大主応力σmax,jの算出対象となる部位ｊとして１つ以上の部位ｊを設定する。
【0164】
そして、応力算出部４４０は、<<回帰係数算出部４２０、回帰係数記憶部４３０、ステップＳ５０２、Ｓ５０３>>の欄で説明した学習時刻ｔ１を、所定の時間隔Δｔ２の推定時刻ｔ２に置き替えた処理を行って、推定時刻ｔ２における部位ｊの状態変数Ξkを算出する（（５２）式、（５７）を参照）。最大主応力σmax,jの算出対象となる部位ｊが複数ある場合、応力算出部４４０は、各推定時刻ｔ２において、各部位ｊの状態変数Ξkをそれぞれ算出する。
【0165】
そして、応力算出部４４０は、各推定時刻ｔ２における部位ｊの状態変数Ξkと、回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jが与えられた（６０）式の線形回帰式、または、回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jが与えられた（６３）式の線形回帰式と、に基づいて、台車枠１６の部位ｊにおける最大主応力σmax,jを、各推定時刻ｔ２において算出する。最大主応力σmax,jの算出対象となる部位ｊが複数ある場合、応力算出部４４０は、各推定時刻ｔ２において、各部位ｊにおける最大主応力σmax,jを算出する。
【0166】
<<出力部４５０、ステップＳ６０３>>
出力部４５０は、応力算出部４４０により算出された台車枠１６の部位ｊにおける最大主応力σmax,jを示す情報を出力する。出力部４５０による情報の出力形態は限定されない。本実施形態の出力部４５０は、第１実施形態の出力部４５０に対し、出力対象の情報が異なるだけであり、本実施形態の出力部４５０が有する機能自体は、第１実施形態の出力部４５０と同じである。
出力部４５０は、応力算出部４４０により算出された台車枠１６の部位ｊにおける最大主応力σmax,jを示す情報を出力する。出力部４５０による情報の出力形態は限定されない。本実施形態の出力部４５０は、第１実施形態の出力部４５０に対し、出力対象の情報が異なるだけであり、本実施形態の出力部４５０が有する機能自体は、第１実施形態の出力部４５０と同じである。
【0167】
＜計算例＞
次に、本実施形態の計算例を説明する。本計算例では、図７における前半の６００個の各時刻において、本実施形態で説明したようにして部位ｊの状態変数Ξk（ξ・p,j (1)、ξp,j (1)、・・・、ξ・p,j (n)、ξp,j (n)）を算出した。これら６００個の各時刻における状態変数Ξkおよび最大主応力を教師データとして、本実施形態で説明したようにして係数行列Ｗj＾を算出した。そして、図７における後半の３００１個の各時刻において、特許文献１に記載の手法と本実施形態の手法とのそれぞれの手法で、台車枠１６の部位ｊの最大主応力を算出した。本計算例では、（６０）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jを成分として含む係数行列Ｗj＾と、（６３）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jを成分として含む係数行列Ｗj＾と、をそれぞれ算出し、それぞれの線形回帰式を用いて台車枠１６の部位ｊの最大主応力をそれぞれ算出した。また、前半の６００個の各時刻においても、前述したようにして算出した係数行列Ｗj＾を用いて本実施形態の手法で台車枠１６の部位ｊの最大主応力を算出した。なお、最大主応力の算出対象の部位ｊは、いずれも同じ部位である。
次に、本実施形態の計算例を説明する。本計算例では、図７における前半の６００個の各時刻において、本実施形態で説明したようにして部位ｊの状態変数Ξk（ξ・p,j (1)、ξp,j (1)、・・・、ξ・p,j (n)、ξp,j (n)）を算出した。これら６００個の各時刻における状態変数Ξkおよび最大主応力を教師データとして、本実施形態で説明したようにして係数行列Ｗj＾を算出した。そして、図７における後半の３００１個の各時刻において、特許文献１に記載の手法と本実施形態の手法とのそれぞれの手法で、台車枠１６の部位ｊの最大主応力を算出した。本計算例では、（６０）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗ2n,jを成分として含む係数行列Ｗj＾と、（６３）式の線形回帰式の回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,jを成分として含む係数行列Ｗj＾と、をそれぞれ算出し、それぞれの線形回帰式を用いて台車枠１６の部位ｊの最大主応力をそれぞれ算出した。また、前半の６００個の各時刻においても、前述したようにして算出した係数行列Ｗj＾を用いて本実施形態の手法で台車枠１６の部位ｊの最大主応力を算出した。なお、最大主応力の算出対象の部位ｊは、いずれも同じ部位である。
【0168】
図９および図１０にその結果を示す。図９は、（６０）式の線形回帰式を用いた結果を示し、図１０は、（６３）式の線形回帰式を用いた結果を示す。図９および図１０において、真値は、特許文献１に記載の手法で算出された台車枠１６の部位ｊの最大主応力である。推定値は、本実施形態の手法で算出された台車枠１６の当該部位ｊの最大主応力である。同時刻における真値および推定値を１つの点としてプロットすることにより、図９（ａ）、図９（ｂ）、図１０（ａ）、および図１０（ｂ）が得られる。
【0169】
図９（ａ）および図１０（ａ）は、前述した前半の６００個の各時刻における、台車枠１６の最大主応力の真値と推定値との関係の一例を示す図である。図９（ｂ）および図１０（ｂ）は、前述した後半の３００１個の各時刻における、台車枠１６の最大主応力の真値と推定値との関係の一例を示す図である。図９（ａ）においては、決定係数Ｒ2は、０．９９７４６３０３３７となり、図９（ｂ）においては、決定係数Ｒ2は、０．９９８９４５９３６６となり、図１０（ａ）においては、決定係数Ｒ2は、０．９９７４０４０１１０となり、図１０（ｂ）においては、決定係数Ｒ2は、０．９９９１４７８４６２となり、本実施形態の手法では、特許文献１に記載の手法と同程度の精度で、台車枠１６の最大主応力を算出することができることが分かる。また、ＣＰＵの処理時間は、特許文献１に記載の手法よりも本実施形態の手法の方が、２２倍以上速く、また、台車枠１６の最大主応力の算出に要する時間は、特許文献１に記載の手法よりも本実施形態の手法の方が、１０００倍以上速くなった。したがって、本実施形態の手法では、特許文献１に記載の手法に比べ、計算時間を大幅に短縮することができることが分かる。
【0170】
＜まとめ＞
以上のように本実施形態では、処理装置４００は、軸箱１７ａ、１７ｂの加速度データを取得する。軸箱１７ａ、１７ｂの加速度データに基づいて、状態変数Ξ（モード座標系における変位ベクトルおよび速度ベクトルの少なくとも一方）を算出する。処理装置４００は、台車枠１６の最大主応力σmax,jを目的変数として含み、状態変数Ξを複数の説明変数として含む線形回帰式に、状態変数Ξを与えることにより、台車枠１６の最大主応力σmax,jを算出する。したがって、第１実施形態と同様に、物体の応力を推定するに際し、計算速度を向上させることができる。
以上のように本実施形態では、処理装置４００は、軸箱１７ａ、１７ｂの加速度データを取得する。軸箱１７ａ、１７ｂの加速度データに基づいて、状態変数Ξ（モード座標系における変位ベクトルおよび速度ベクトルの少なくとも一方）を算出する。処理装置４００は、台車枠１６の最大主応力σmax,jを目的変数として含み、状態変数Ξを複数の説明変数として含む線形回帰式に、状態変数Ξを与えることにより、台車枠１６の最大主応力σmax,jを算出する。したがって、第１実施形態と同様に、物体の応力を推定するに際し、計算速度を向上させることができる。
【0171】
＜変形例＞
本実施形態では、台車枠１６の各メッシュの節点の位置における応力（本実施形態では最大主応力）を各学習時刻ｔ１において算出し、学習時刻ｔ１における部位ｊの教師データｐとする場合を例示した。
しかしながら、教師データに含める目的変数は、応力の計算値ではなく、測定値であってもよい（以下の説明では、応力の測定値を、必要に応じて応力測定値と称する）。
本実施形態では、台車枠１６の各メッシュの節点の位置における応力（本実施形態では最大主応力）を各学習時刻ｔ１において算出し、学習時刻ｔ１における部位ｊの教師データｐとする場合を例示した。
しかしながら、教師データに含める目的変数は、応力の計算値ではなく、測定値であってもよい（以下の説明では、応力の測定値を、必要に応じて応力測定値と称する）。
【0172】
図１１は、台車枠１６の部位ｊの応力測定値と時間との関係の一例を示す図である。図１１では、台車枠１６の部位ｊが横梁と心皿とのつなぎ部である場合を例示する。また、ここでは、横梁と心皿とのつなぎ部に貼付した歪みゲージを用いて、当該つなぎ部の応力測定値を得た。また、ここでは、教師データに含める目的変数が応力測定値であることから、以下の（６６）式に示すように、（６３）式の左辺に示すσmax,jを、σobs,jに置き換えて表記する（（６６）式の右辺は（６３）式の右辺と同じである）。
【0173】
図１１における前半の１０００個の各時刻ｋの応力測定値のデータと、当該時刻のモード座標系における変位ベクトル[ξ]と、を教師データとして用いて、本実施形態で説明した手法で（６６）式の係数行列Ｗj＾（回帰係数ｗ0,j、ｗ1,j、・・・、ｗn,j）を決定した。そして、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊにおける応力σobs,jを、学習後の（６６）式を用いて計算（推定）した。図１１における後半の３００１個の各時刻における計算値（推定値）と、当該各時刻における測定値との関係を図１２に示す。ここでは、モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］、外力ベクトル｛ｆξ｝を（５１）式に与えて、（５１）式を解くことにより、モード座標系における台車枠１６の部位ｊの変位ベクトル｛ξ｝を算出する場合を例示する。
【0174】
【数35】
【0175】
図１２において、推定値は、学習後の（６６）式を用いて算出した、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力σobs,jである。真値は、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力測定値である。同時刻における測定値および計算値を１つの点としてプロットすることにより、図１２が得られる。図１２に示すように、決定係数Ｒ2は、０．５９９となり、教師データとして応力測定値を使っても良好な推定精度が得られることが分かる。このようにする場合、図５のステップＳ５１０において、データ取得部４１０は、台車枠１６の応力測定値をさらに取得する。
【0176】
ここでは、（６３）式の線形回帰式を例示して、教師データに含める目的変数が応力測定値でもよいことを説明した。（６０）式においても、教師データに含める目的変数を応力測定値としてもよい。
【0177】
また、特許文献２に記載のように、モード座標系における台車枠１６の変位ベクトル｛ξ｝の推定精度を向上させるために、データ同化を行うフィルタを用いてもよい。すなわち、応力算出部４４０は、物理座標系における台車枠１６の変位ベクトル[ｕ]を観測変数ｙとすると共に、モード座標系における台車枠１６の変位ベクトル[ξ]および速度ベクトル[ξ・]を状態変数Ξとし、観測変数ｙの測定値と（５１）式に基づく計算値との差異が小さく（最小に）なるように、未観測の変数（状態変数Ξ）の推定値を、データ同化を行うフィルタを用いた演算を行うことにより算出してもよい。観測変数ｙの測定値（物理座標系における台車枠１６の変位ベクトル[ｕ]）は、加速度センサ２２ａ、２２ｂで測定された台車枠１６の加速度データに基づいて算出される。
【0178】
観測変数ｙ∈Ｒ3mは、例えば、以下の（５７）式のように表される。ｍ∈Ｚは、観測位置の数である。また、観測方程式は、以下の（６８）式で表される。
【0179】
【数36】
【0180】
図２に示す例では、加速度センサ２２ａ、２２ｂの数は２であるので、観測位置の数ｍは２である。～は、測定量であることを表す。尚、測定量には、直接測定される量だけではなく、他の測定量を用いて導出される量も含まれるものとする。
（６７）式において、ｕ～１，ｘ、ｕ～１，ｙ、ｕ～１，ｚ、・・・ｕ～ｍ，ｘ、ｕ～ｍ，ｙ、ｕ～ｍ，ｚのカンマ（，）の前の１、・・・、ｍの値をとる変数は、観測位置（図２に示す例では、加速度センサ２２ａ、２２ｂの取り付け位置）を表す。また、ｕ～１，ｘ、ｕ～１，ｙ、ｕ～１，ｚ、・・・ｕ～ｍ，ｘ、ｕ～ｍ，ｙ、ｕ～ｍ，ｚのカンマ（，）の後のｘ、ｙ、ｚは、それぞれ、ｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分であることを表す。従って、例えば、ｕ～１，ｘは、変数＝１の観測位置に対応する加速度センサの取り付け位置における変位のｘ軸方向成分の測定量を表す。
なお、ｕ～は、（６７）式においてｕの上に～が付されているものに対応する。また、（６８）式において、Ｖは、観測ノイズである。
（６７）式において、ｕ～１，ｘ、ｕ～１，ｙ、ｕ～１，ｚ、・・・ｕ～ｍ，ｘ、ｕ～ｍ，ｙ、ｕ～ｍ，ｚのカンマ（，）の前の１、・・・、ｍの値をとる変数は、観測位置（図２に示す例では、加速度センサ２２ａ、２２ｂの取り付け位置）を表す。また、ｕ～１，ｘ、ｕ～１，ｙ、ｕ～１，ｚ、・・・ｕ～ｍ，ｘ、ｕ～ｍ，ｙ、ｕ～ｍ，ｚのカンマ（，）の後のｘ、ｙ、ｚは、それぞれ、ｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分であることを表す。従って、例えば、ｕ～１，ｘは、変数＝１の観測位置に対応する加速度センサの取り付け位置における変位のｘ軸方向成分の測定量を表す。
なお、ｕ～は、（６７）式においてｕの上に～が付されているものに対応する。また、（６８）式において、Ｖは、観測ノイズである。
【0181】
（６８）式に示す観測方程式により、観測変数ｙは、状態変数Ξと関連付けられる。観測行列Ｈ∈Ｒ3m×2nは、（４８）式より、以下の（６９）式で表される。
【0182】
【数37】
【0183】
（６９）式において、φ1,x
(1)、φ1,x
(2)、・・・、φ1,x
(n)、φ1,y
(1)、φ1,y
(2)、・・・、φ1,y
(n)、φ1,z
(1)、φ1,z
(2)、・・・、φ1,z
(n)、・・・、φm,x
(1)、φm,x
(2)、・・・、φm,x（ｎ）、φm,y
(1)、φm,y
(2)、・・・、φm,y
(n)、φm,z
(1)、φm,z
(2)、・・・、φm,z
(n)のカンマ（，）の前の１、・・・、ｍの値をとる変数は、観測位置（図２に示す例では、加速度センサ２２ａ、２２ｂの取り付け位置）を表す。φ1,x
(1)、φ1,x
(2)、・・・、φ1,x
(n)、φ1,y
(1)、φ1,y
(2)、・・・、φ1,y
(n)、φ1,z
(1)、φ1,z
(2)、・・・、φ1,z
(n)、・・・、φm,x
(1)、φm,x
(2)、・・・、φm,x（ｎ）、φm,y
(1)、φm,y
(2)、・・・、φm,y
(n)、φm,z
(1)、φm,z
(2)、・・・、φm,z
(n)のカンマ（，）の後のｘ、ｙ、ｚは、それぞれ、ｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分であることを表す。従って、例えば、φ1,x
(1)は、１次の固有振動モード、変数＝１の観測位置に対応する加速度センサの取り付け位置、およびｘ軸方向成分に対応するモード行列［φ］の成分である。
【0184】
（５２）式、（６８）式、および（６９）式に示すように、状態変数Ξのうち、モード座標系における台車枠１６の変位ξ(1)、・・・、ξ(n)のみに、モード行列［φ］の成分が乗算される。具体的に、モード座標系における台車枠１６の変位ξ(1)、・・・、ξ(n)に乗算されるモード行列［φ］の成分は、当該モード座標系における台車枠１６の変位と固有振動モードの次数が同じモード行列［φ］の成分であり、且つ、当該モード座標系における台車枠１６の変位が生じる観測位置と同じ観測位置に対応するモード行列［φ］の成分である。
【0185】
（５３）式に示す状態方程式（ただし、（５３）式の右辺にシステムノイズが加算される）と、（６７）式に示す観測方程式と、加速度センサ２２ａ、２２ｂで測定された加速度のデータから導出された観測変数の測定値と、を用いて、データ同化を行うフィルタに基づく演算を行うことにより、モード座標系における台車枠１６の変位ベクトル｛ξ｝（＝［ξ(1)、・・・、ξ(n)］T）の推定値（モード座標系における台車枠１６の変位の分布の推定値）が算出される。この場合、台車枠１６の加速度データが、物体の振動に応じて変化する第１物理量の測定データの一例になる。
【0186】
データ同化を行うフィルタの一例として、線形カルマンフィルタがある。線形カルマンフィルタでは、観測変数の測定値と計算値との誤差または当該誤差の期待値が小さく（最小に）なるカルマンゲインを求め、そのときの未観測の変数（状態変数）の値を求める。線形カルマンフィルタ自体は、公知の技術で実現することができる。なお、特許文献２に記載されているように、観測変数の測定値と計算値との誤差が最小または当該誤差の期待値が最小になるように状態変数の推定値を導出するフィルタ（即ち、データ同化を行うフィルタ）として線形カルマンフィルタ以外のフィルタ（例えば、粒子フィルタ）を用いてもよい。なお、観測変数の測定値と計算値との誤差としては、例えば、観測変数の測定値と計算値との二乗誤差が挙げられる。データ同化を行うフィルタを用いる手法自体の詳細は、特許文献２に記載されているので、ここでは、その詳細な説明を省略する。
【0187】
図１３（ａ）、図１４（ａ）、図１５（ａ）、図１６（ａ）、図１７（ａ）、図１８（ａ）は、各々、台車枠１６のバネ帽１位、バネ帽２位、バネ帽３位、バネ帽４位、モータ座１軸、モータ座２軸の上下方向（ｚ軸方向）の変位の測定値１３１０、１４１０、１５１０、１６１０、１７１０、１８１０と第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０との関係の一例を示す図である。第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０は、モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］、外力ベクトル｛ｆξ｝を（５１）式に与えて、（５１）式を解くことにより算出された、モード座標系における台車枠１６の変位ベクトル｛ξ｝を物理座標系の変位ベクトル[ｕ]に変換することにより得られた計算値である。
【0188】
不図示のバネ帽は、軸バネ１９ａ、１９ｂを収容する。バネ帽１位は、軸バネ１９ａに対して設けられるバネ帽であり、バネ帽２位は、車軸１５ａの両側のうち、軸バネ１９ａが配置されている側とは反対側の不図示の軸バネに対して設けられるバネ帽である。バネ帽３位は、、軸バネ１９ｂに対して設けられるバネ帽であり、バネ帽２位は、車軸１５ｂの両側のうち、軸バネ１９ｂが配置されている側とは反対側の不図示の軸バネに対して設けられるバネ帽である。
【0189】
不図示のモータ座（モータブラケット）は、台車枠１６の不図示の横梁に設けられ、モータを固定するためのものである。モータ座１軸は、輪軸１７ａ側の横梁に対して設けられるモータ座であり、モータ座２軸は、輪軸１７ｂ側の横梁に対して設けられるモータ座である。
【0190】
図１３（ｂ）、図１４（ｂ）、図１５（ｂ）、図１６（ｂ）、図１７（ｂ）、図１８（ｂ）は、各々、台車枠１６のバネ帽１位、バネ帽２位、バネ帽３位、バネ帽４位、モータ座１軸、モータ座２軸の上下方向（ｚ軸方向）の変位の測定値１３１０、１４１０、１５１０、１６１０、１７１０、１８１０と第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０との関係の一例を示す図である。第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０は、前述したように、特許文献２に記載のようにデータ同化を行うフィルタを用いて算出されたード座標系における台車枠１６の変位ベクトル[ξ]を物理座標系における変位ベクトル[ｕ]に変換することにより得られた計算値である。
【0191】
図１３～図１８に示すように、第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０よりも、第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０の方が、測定値１３１０、１４１０、１５１０、１６１０、１７１０、１８１０に近く、各部位の変位を高精度に計算できることが分かる。尚、図１３～図１８において、測定値１３１０、１４１０、１５１０、１６１０、１７１０、１８１０をグレー（薄い濃度）で示し、第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０および第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０を黒（濃い濃度）で示す。このような表記により、測定値１３１０、１４１０、１５１０、１６１０、１７１０、１８１０のうち、第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０および第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０と同じ値の部分が第１計算値１３２０、１４２０、１５２０、１６２０、１７２０、１８２０および第２計算値１３３０、１４３０、１５３０、１６３０、１７３０、１８３０に隠れて見えなくなっている。
【0192】
データ同化を行うフィルタを用いることにより、モード座標系における変位ベクトル[ξ]は、データ同化を行うフィルタのアルゴリズムに基づいて更新され、各時刻での座標系における変位ベクトル[ξ]として算出される。図１１における前半の１０００個の各時刻においてこのようにして算出された台車枠１６のモード座標系における変位ベクトル[ξ]と、図１１に示す当該各時刻における応力測定値と、を教師データとして用いて、本実施形態で説明した手法で（６６）式の係数行列Ｗj＾（回帰係数（ｗ0,j、ｗ1,j、・・・、ｗn,j）を決定した。
【0193】
図１１における後半の３００１個の各時刻における計算値（推定値）と、当該各時刻における測定値との関係を図１９に示す。
図１９において、推定値は、学習後の（６６）式を用いて算出した、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力σobs,jである。真値は、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力測定値である。同時刻における測定値および計算値を１つの点としてプロットすることにより、図１９が得られる。図１９に示すように、決定係数Ｒ2は、０．６８７となり、図１２に示す結果（決定係数Ｒ2＝０．６５６）に比べて、データ同化を行うフィルタを用いて決定したモード座標系における変位ベクトル[ξ]を説明変数として用いることにより、台車枠１６の応力の推定精度がより向上する結果が得られる。
図１９において、推定値は、学習後の（６６）式を用いて算出した、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力σobs,jである。真値は、図１１における後半の３００１個の各時刻における台車枠１６の部位ｊの応力測定値である。同時刻における測定値および計算値を１つの点としてプロットすることにより、図１９が得られる。図１９に示すように、決定係数Ｒ2は、０．６８７となり、図１２に示す結果（決定係数Ｒ2＝０．６５６）に比べて、データ同化を行うフィルタを用いて決定したモード座標系における変位ベクトル[ξ]を説明変数として用いることにより、台車枠１６の応力の推定精度がより向上する結果が得られる。
【0194】
ここでは、教師データに含まれる目的変数が応力測定値である場合を例示して、データ同化を行うフィルタを用いて算出されたード座標系における台車枠１６の変位ベクトル[ξ]を説明変数として用いる場合を例示した。しかしながら、教師データに含まれる目的変数は、本実施形態で説明したように計算値であってもよい。また、（６０）式の線形回帰式においても同様の結果を得ることができる。したがって、（６０）式の線形回帰式において、データ同化を行うフィルタを用いて算出されたード座標系における台車枠１６の変位ベクトル[ξ]を説明変数として用いてもよい。
【0195】
また、第１実施形態では、物理座標系における座標（ｘ，ｙ，ｚ）をアフィン変換する場合を例示した。しかしながら、第１実施形態において、物理座標系ではなくモード座標系においてアフィン変換を行ってもよい。このようにする場合、（３）式を、例えば、以下の（７０）式に置き替えて、第１実施形態で説明した台車枠１６の部位ｊにおける最大主応力σmax,jを表す線形回帰式を求めてもよい。
【0196】
【数38】
【0197】
（７０）式において、φm,x
(r)、φm,y
(r)、φm,z
(r)、ξm,t+1
(r)、ξm,t
(r)のカンマ（，）の前のｍは、観測位置（第１実施形態の部位ｉに対応する位置）を表す。φm,x
(r)、φm,y
(r)、φm,z
(r)のカンマ（，）の後のｘ、ｙ、ｚは、それぞれ、ｘ軸方向成分、ｙ軸方向成分、ｚ軸方向成分であることを表す。ξm,t+1
(r)、ξm,t
(r)のカンマ（，）の後のｔ＋１、ｔは、時刻ｔ＋１、ｔを表す。時刻ｔ＋１は、時刻ｔに時間ステップΔｔを加算した次の時刻である。モード座標系における台車枠１６の変位ベクトル｛ξ｝の推定精度を向上させるために、本変形例で示したようにデータ同化を行うフィルタを用いて算出するのが好ましい。しかしながら、モード座標系における台車枠１６の質量行列［Ｍξ］、粘性行列［Ｃξ］、剛性行列［Ｋξ］、外力ベクトル｛ｆξ｝を（５１）式に与えて、（５１）式を解くことにより、モード座標系における台車枠１６の部位ｊの変位ベクトル｛ξ｝を算出してもよい。
【0198】
なお、本実施形態においても、第１実施形態で説明した種々の変形例のうち、アフィン変換行列に特化した変形例を除く変形例を採用してもよい。
【0199】
（その他の実施形態）
なお、以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体及び前記プログラム等のコンピュータプログラムプロダクトも本発明の実施形態として適用することができる。記録媒体としては、例えば、フレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ－ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。また、本発明の実施形態は、ＰＬＣ（Programmable Logic Controller）により実現されてもよいし、ＡＳＩＣ（Application Specific Integrated Circuit）等の専用のハードウェアにより実現されてもよい。
また、以上説明した本発明の実施形態は、いずれも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
なお、以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体及び前記プログラム等のコンピュータプログラムプロダクトも本発明の実施形態として適用することができる。記録媒体としては、例えば、フレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ－ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。また、本発明の実施形態は、ＰＬＣ（Programmable Logic Controller）により実現されてもよいし、ＡＳＩＣ（Application Specific Integrated Circuit）等の専用のハードウェアにより実現されてもよい。
また、以上説明した本発明の実施形態は、いずれも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
【0200】
なお、以上の実施形態の開示は、例えば以下のようになる。
（開示１）
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得手段と、
前記データ取得手段により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出手段と、
を有し、
前記応力算出手段は、前記データ取得手段により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理装置。
（開示２）
前記複数の説明変数は、変換前の第２物理量と、当該変換前の第２物理量に乗算される係数と、に基づいて、変換後の前記第２物理量を算出する計算式における、前記係数、または、前記変換後の第２物理量を含み、
前記第２物理量は、位置、変位、および速度のうち、少なくとも１つを含む、開示１に記載の処理装置。
（開示３）
前記計算式は、前記物体の変位前の物理座標系の位置の座標と、当該座標に乗算される変換行列と、に基づいて、前記物体の変位後の物理座標系の位置を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、開示２に記載の処理装置。
（開示４）
前記計算式は、物理座標系における前記物体の振動を記述する運動方程式を固有値解析することにより算出される固有ベクトルと、当該固有ベクトルに乗算されるモード座標系における前記物体の変位と、に基づいて、物理座標系における前記物体の変位を算出する計算式と、前記固有ベクトルと、モード座標系における前記物体の速度と、に基づいて、物理座標系における前記物体の速度を算出する計算式と、を含み、
前記複数の説明変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、開示２に記載の処理装置。
（開示５）
前記計算式は、前記物体のモード座標系における前記物体の変位と、当該座標に乗算される変換行列と、に基づいて、当該変位が生じている時刻の次の時刻での、前記物体のモード座標系における前記物体の変位を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、開示２に記載の処理装置。
（開示６）
前記応力算出手段は、モード座標系における前記物体の振動を表す運動方程式に基づいて構成される状態方程式と、観測変数と状態変数との関係を表す観測方程式と、前記第１物理量の測定データから算出される前記観測変数の測定値と、に基づいて、データ同化を行うフィルタを用いた演算を行うことにより、前記状態変数の推定値を決定し、
前記状態変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、開示４または５に記載の処理装置。
（開示７）
前記物体は、鉄道車両の台車枠を含む、開示１～６のいずれか１つに記載の処理装置。
（開示８）
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得工程と、
前記データ取得工程により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出工程と、
を有し、
前記応力算出工程は、前記データ取得工程により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理方法。
（開示９）
開示１～７のいずれか１つに記載の処理装置の各手段としてコンピュータを機能させるためのプログラム。
（開示１）
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得手段と、
前記データ取得手段により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出手段と、
を有し、
前記応力算出手段は、前記データ取得手段により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理装置。
（開示２）
前記複数の説明変数は、変換前の第２物理量と、当該変換前の第２物理量に乗算される係数と、に基づいて、変換後の前記第２物理量を算出する計算式における、前記係数、または、前記変換後の第２物理量を含み、
前記第２物理量は、位置、変位、および速度のうち、少なくとも１つを含む、開示１に記載の処理装置。
（開示３）
前記計算式は、前記物体の変位前の物理座標系の位置の座標と、当該座標に乗算される変換行列と、に基づいて、前記物体の変位後の物理座標系の位置を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、開示２に記載の処理装置。
（開示４）
前記計算式は、物理座標系における前記物体の振動を記述する運動方程式を固有値解析することにより算出される固有ベクトルと、当該固有ベクトルに乗算されるモード座標系における前記物体の変位と、に基づいて、物理座標系における前記物体の変位を算出する計算式と、前記固有ベクトルと、モード座標系における前記物体の速度と、に基づいて、物理座標系における前記物体の速度を算出する計算式と、を含み、
前記複数の説明変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、開示２に記載の処理装置。
（開示５）
前記計算式は、前記物体のモード座標系における前記物体の変位と、当該座標に乗算される変換行列と、に基づいて、当該変位が生じている時刻の次の時刻での、前記物体のモード座標系における前記物体の変位を算出する計算式を含み、
前記複数の説明変数は、前記計算式における変換行列の成分のうち、線形変換を行うための成分を含む、開示２に記載の処理装置。
（開示６）
前記応力算出手段は、モード座標系における前記物体の振動を表す運動方程式に基づいて構成される状態方程式と、観測変数と状態変数との関係を表す観測方程式と、前記第１物理量の測定データから算出される前記観測変数の測定値と、に基づいて、データ同化を行うフィルタを用いた演算を行うことにより、前記状態変数の推定値を決定し、
前記状態変数は、モード座標系における前記物体の変位および速度の少なくとも一方を含む、開示４または５に記載の処理装置。
（開示７）
前記物体は、鉄道車両の台車枠を含む、開示１～６のいずれか１つに記載の処理装置。
（開示８）
物体の振動に応じて変化する第１物理量の測定データを含むデータを取得するデータ取得工程と、
前記データ取得工程により取得されたデータと、前記物体の応力を目的変数として含み、前記物体の応力に対する影響因子を複数の説明変数として含む線形回帰式と、に基づいて、前記物体の応力を算出する応力算出工程と、
を有し、
前記応力算出工程は、前記データ取得工程により取得された前記第１物理量の測定データに基づいて、前記複数の説明変数を算出する、処理方法。
（開示９）
開示１～７のいずれか１つに記載の処理装置の各手段としてコンピュータを機能させるためのプログラム。
【符号の説明】
【0201】
１１ 車体
１２ａ、１２ｂ 台車
１３ａ～１３ｄ 輪軸
１４ａ～１４ｄ 車輪
１５ａ～１５ｄ 車軸
１６ 台車枠
１７ａ、１７ｂ 軸箱
１８ａ、１８ｂ モノリンク
１９ａ、１９ｂ 軸バネ
２１ａ、２１ｂ 加速度センサ
２２ａ、２２ｂ 加速度センサ
３１ａ、３１ｂ 着力点
３２ａ、３２ｂ 着力点
４００ 処理装置
４１０ データ取得部
４２０ 回帰係数算出部
４３０ 回帰係数記憶部
４４０ 応力算出部
４５０ 出力部
１２ａ、１２ｂ 台車
１３ａ～１３ｄ 輪軸
１４ａ～１４ｄ 車輪
１５ａ～１５ｄ 車軸
１６ 台車枠
１７ａ、１７ｂ 軸箱
１８ａ、１８ｂ モノリンク
１９ａ、１９ｂ 軸バネ
２１ａ、２１ｂ 加速度センサ
２２ａ、２２ｂ 加速度センサ
３１ａ、３１ｂ 着力点
３２ａ、３２ｂ 着力点
４００ 処理装置
４１０ データ取得部
４２０ 回帰係数算出部
４３０ 回帰係数記憶部
４４０ 応力算出部
４５０ 出力部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】