JP6905296B1 - 制御対象の動作を模擬する計算装置、計算装置において状態方程式を解く演算の実行方法、及びプログラム - Google Patents

Abstract

【課題】制御対象の動作を数学的にモデル化した状態方程式を解く演算を実行することにより制御対象の動作を模擬し、それにより制御装置の設計及び検証を支援する計算装置の技術に関し、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とする。【解決手段】固有値計算手段３０２は、変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する。固有値変更手段１２０１は、固有値包含対角行列内の各固有値が所定の条件を満たす場合に、固有値の実部又は虚数部の一方又は両方を所定の定数に変更し、修正固有値包含対角行列を出力する。方程式変換手段である行列算出手段１２０２は、線形状態方程式を修正固有値包含対角行列に基づく形式に変換し、修正線形状態方程式を出力する。状態方程式計算手段１２０３は、修正線形状態方程式を差分方程式の形式に変換し、差分方程式に対する解法演算を実行することにより、制御装置への出力信号を生成する。【選択図】図１２

Description

本発明は、制御対象の動作を数学的にモデル化した状態方程式を解く演算を実行することにより制御対象の動作を模擬し、それにより制御装置の設計及び検証を支援する計算装置、そのような計算装置において状態方程式を解く演算の実行方法、及びプログラムに関する。

制御対象を接続することにより、その制御対象に対してフィードバック制御を行う制御装置が知られている。その制御対象は、機械装置である場合もあり、或いは電子回路である場合もある。例えば、自動車には、ＥＣＵ（Ｅｌｅｃｔｒｏｎｉｃｓ Ｃｏｎｔｒｏｌ Ｕｎｉｔ：電子制御ユニット）と呼ばれる制御装置が、数十個搭載されている。夫々のＥＣＵは、例えばエンジン、トランスミッション、又はサスペンションなどの機械装置である制御対象に対して、フィードバック制御を行う。また、ＥＣＵは、例えばダッシュボードパネルに各種情報を表示する液晶計器類をはじめとする電子回路である制御対象に対して、フィードバック制御を行う。これらの制御装置は、制御対象を接続して初めて機能するという特性を持つ。従って、制御装置の開発においては、かつては、実際に制御対象を接続して実験／検証を行う必要があった。しかし、開発段階の制御装置を実際の制御対象に接続して検証を行った場合、制御装置による制御が不完全であることにより、制御対象が破壊されてしまうという事故が起こり得た。また、例えば月着陸船の制御装置のように、実際の使用環境での実験／検証が不可能な制御装置の開発においては、制御対象を実際に接続してその制御装置を開発することができなかった。或いは、開発段階において制御対象の実機が存在しない場合においても、制御対象を実際に接続してその制御装置を開発することができなかった。

このような問題を解決するために、従来、ＨＩＬＳ（Ｈａｒｄｗａｒｅ Ｉｎ ｔｈｅ Ｌｏｏｐ Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置が知られている（例えば特許文献１）。ＨＩＬＳでは、機械装置又は電子回路である制御対象が計算装置におけるソフトウェア処理によって模擬される。この場合、制御対象の動作は、例えば常微分方程式である状態方程式によって構成される数学モデルで表現される。また、計算装置は、電圧又は電流等の実際の信号を入出力する物理的入出力手段を装備する。そして、計算装置のこれらの物理的入出力手段に、制御装置の実際の物理的入出力手段が接続される。計算装置は、状態方程式を解く演算を実時間で実行することにより制御対象の動作を模擬する。そして、計算装置は実際の信号を制御装置に対して実時間で入出力することにより、ハードウェアである制御装置を実際に動作させてテストする。即ち、ＨＩＬＳは、リアルタイムシミュレータとして動作する。ＨＩＬＳにより、制御対象の実機を用意することなく、ＥＣＵなどの制御装置の設計及び検証を支援することが可能となる。

更に、ＳＩＬＳ（Ｓｏｆｔｗａｒｅ Ｉｎ ｔｈｅ Ｌｏｏｐ Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置も知られている（例えば特許文献１）。ＳＩＬＳでは、制御装置自体、及びその物理的入出力手段も計算装置によって模擬される。そして、ＳＩＬＳは、制御装置と制御対象の全体の動作を、計算装置が実行するソフトウェア処理によって模擬する。

ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、一般的に常微分方程式である状態方程式として表現することができる。特に、変数ベクトルに第１の係数行列を乗算し、入力ベクトルに第２の係数行列を乗算した各乗算結果の和が、変数ベクトルの１階の導関数に等しい線形状態方程式は、制御対象の挙動を表す数学モデルとして有用であり、非線形の場合であっても上記形式に近似して用いられることが多い。ここで、「変数ベクトル」とは、常微分方程式の解をいう。

ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、線形状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実時間で入出力する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。

連続系の状態方程式を解くための計算方法としては、その状態方程式を離散信号形式で表現することにより差分方程式として計算する数値計算手法が良く知られている。この解法の代表例としては、Ｒｕｎｇｅ−Ｋｕｔｔａ法、及びＴａｙｌｏｒ展開法などが知られている。

ところで、線形状態方程式における係数行列の要素が或る一定時間において不変な定数であれば、その方程式は、所謂ＬＴＩ（Ｌｉｎｅａｒ Ｔｉｍｅ Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルは、ＬＴＩシステムを仮定しても十分であることが多い。このようなＬＴＩシステムの状態方程式においては、係数行列を対角要素が互いに異なる固有値である対角行列に変換することができる。

このような対角行列に含まれる固有値は、一般に複素数である。そして、複数の固有値の大きさが著しく異なる状態方程式は、解が周波数の比較的高い成分と低い成分を含む、いわゆるＳｔｉｆｆ状態方程式と呼ばれる（例えば特許文献２）。

特許第５６９２７３９号公報 特開２００５−０２５６５１号公報

上述のようなＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、前述したＲｕｎｇｅ−Ｋｕｔｔａ法又はＴａｙｌｏｒ展開法などによる差分方程式に変換して演算する場合について考える。このケースでは、差分方程式における離散時間の計算ステップと対角行列の固有値が適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（−∞）に発散することが知られている。特に、前述したＳｔｉｆｆ状態方程式を解く場合に、絶対値の大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくしなければならない。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に対しては、その計算ステップは過剰に小さくなりすぎる。この結果、計算ステップを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。

このような問題を解決するために、更に従来、複素平面のいわゆる左半平面全体を計算安定領域とする、後退Ｅｕｌｅｒ法又は台形法（Ｔｕｓｔｉｎ法）と呼ばれる解法が知られている。しかしながら、これらの解法ではＬＴＩシステムの場合であっても逆行列を計算する必要がある。一般に逆行列の計算負荷は、行列が“ｎ×ｎ”要素のサイズを有する場合、ｎの３乗に比例すると言われている。従って、ｎが大きなシステムの場合には、計算負荷が極めて高くなる。このため、Ｅｕｌｅｒ法や台形法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。

更に、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次、台形法は２次の近似精度となる。このため、固有値λの絶対値が大きい場合は、これらの解法は、安定動作範囲が広いため発散しないで計算できるものの、計算精度が著しく落ちるというあまり知られていない事実がある。

以上の考察からわかるように、計算装置において絶対値が大きな固有値を有する方程式を適切に解く演算を実行することは、従来困難であった。

そこで、本発明は、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とすることを目的とする。

態様の一例は、制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、制御対象の動作を模擬する計算装置であって、制御装置への出力を表す変数ベクトル項Ｙと制御装置からの入力を表す入力ベクトル項Ｕを含み、変数ベクトル項Ｘ及び入力ベクトル項Ｕに夫々乗算される各係数行列Ａ，Ｂの要素が一定時間において不変な定数とみなせる線形時間不変システムに制御対象をモデル化した線形状態方程式を、状態方程式として設定する方程式設定手段と、変数ベクトル項Ｘの係数行列Ａを、対角要素が固有値である対角行列Λ又はジョルダン標準形行列Ｊの何れかの行列である固有値包含対角行列に変換する固有値計算手段と、固有値包含対角行列内の各固有値が所定の条件を満たす場合に、固有値の実部又は虚数部の一方又は両方を所定の定数に変更することにより、固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更手段と、線形状態方程式を修正固有値包含対角行列に基づく形式に変換することにより、線形状態方程式に対応する修正線形状態方程式を出力する方程式変換手段と、修正線形状態方程式を差分方程式の形式に変換し、差分方程式に対する解法演算を実行することにより、制御装置への出力信号を生成する状態方程式計算手段と、を備える。

本発明によれば、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、ＨＩＬＳまたはＳＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。

Ｒｕｎｇｅ−Ｋｕｔｔａ法において安定に計算できる固有値λと計算ステップhの積ｈλの複素平面領域を説明する図である。 Ｔａｙｌｏｒ展開法において安定に計算できる固有値λと計算ステップｈの積ｈλの領域を説明する図である。 本発明の第１の実施形態であるＨＩＬＳ計算装置の実施形態を示すブロック図である。 制御対象の機構例を示す図である。 制御対象のブロック線図の例を示す図である。 固有値実部変更手段及び固有値虚数部変更手段の動作例を示すフローチャートである。 本発明の第２の実施形態であるＳＩＬＳ計算装置の実施形態を示すブロック図である。 固有値λ＝±πｉを異なる計算手法で計算した時系列結果を説明する図である。 固有値λ＝±π／２ｉを異なる計算手法で計算した時系列結果を説明する図である。 指数関数の時間応答を異なる計算手法で計算した時系列結果を説明する図である。 複素数の固有値を異なる計算手法で計算した時系列結果を説明する図である。 本発明の第３の実施形態であるＨＩＬＳ計算装置の実施形態を示すブロック図である。 第３の実施形態における固有値λの変更方法を説明する図である。 固有値変更手段の動作例を示すフローチャートである。 ４次のＴａｙｌｏｒ展開法におけるゼロ点の複素平面上での位置を説明する図である。 固有値が−０．２７±２．５０ｉの場合において４次のＴａｙｌｏｒ展開法で計算した時の時系列計算結果と、第４の実施形態による改善例を説明する図である。

以下、本発明を実施するための形態（以下「本実施形態」と記載）について図面を参照しながら詳細に説明する。以下の説明においては、本実施形態について詳細に説明する前に、まず状態方程式に関する一般的な技術背景について説明する。そして、その一般技術の問題点を具体的に明らかにした上で、本実施形態について詳細に説明する。

前述したように、ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、任意数の要素を持つ状態変数ベクトルＸ、時間を表す変数ｔ、及びベクトル関数Ｆ（ｔ，Ｘ）を用いて、下記（１）式の常微分方程式として表現される。

上記（１）式の常微分方程式は、公知の技術である動作点による区分線形化手法などにより、Ａ、Ｂ行列と、外力を表す入力ベクトルＵ（ｔ）とを用いた、線形状態方程式と呼ばれる下記（２）式の形式で表現することができる。

このように、線形状態方程式は、変数ベクトルＸ（ｔ）に係数行列Ａを乗算し、入力ベクトルＵ（ｔ）に係数ベクトルＢを乗算した各乗算結果ＡＸ（ｔ）とＢＵ（ｔ）の和が、変数ベクトルの１階の導関数

に等しい、という形を有する。

ＨＩＬＳまたはＳＩＬＳにおいて（２）式の形式の状態方程式を解く演算を効率的に実行することは重要である。特に、ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、（２）式の状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実現時間で供給する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。

上述の（２）式に示されるようなＡ、Ｂ行列で表現される連続系方程式を解くための計算方法としては、（２）式の形式の状態方程式を離散信号形式で表現することにより差分方程式として計算する数値計算手法が知られている。前述したように、この解法の代表例として従来、Ｒｕｎｇｅ−Ｋｕｔｔａ法、又はＴａｙｌｏｒ展開法などが知られている。

ところで、上記（２）式の各係数行列Ａ、Ｂは必ずしも定数行列である必要はない。しかしながら、各係数行列Ａ、Ｂの要素が一定時間において不変な定数であれば、（２）式は、所謂ＬＴＩ（Ｌｉｎｅａｒ Ｔｉｍｅ Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの線形状態方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルはＬＴＩシステムを仮定しても十分であることが多い。係数行列Ａ、Ｂが定数であるＬＴＩシステムでは、（２）式は下記（３）、（４）式に示されるように固有値分解することができる。

上記（４）式に示されるように、Λは対角要素が固有値である対角行列である。従って、（２）式は、スカラ連立方程式として扱うことができる。

上記（４）式の固有値λ_ｓ（１≦ｓ≦ｎ）は、一般に複素数である。そして、前述したように、複数の固有値λ_ｓの大きさが著しく異なる状態方程式は、いわゆるＳｔｉｆｆ状態方程式と呼ばれる。

なお、係数行列Ａを上記（４）式のように対角行列に分解ができない場合もある。その場合には、後述するように、係数行列Ａはジョルダン（Ｊｏｒｄａｎ）標準形行列に分解することができる。

上述の（３）、（４）式に示されるＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、前述したＲｕｎｇｅ−Ｋｕｔｔａ法又はＴａｙｌｏｒ展開法などの差分方程式に変換して演算する場合について考える。このケースでは、差分方程式における離散時間の計算ステップｈと、Ａ行列の固有値λの値とが適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（−∞）に発散することが知られている。逆に発散せず安定に状態方程式を解くことができる条件は、計算アルゴリズムごとに異なる。

図１は、Ｒｕｎｇｅ−Ｋｕｔｔａ法において安定に計算できる固有値λと計算ステップｈの積ｈλの複素平面領域を説明する図である。図１において、横軸は積ｈλの値の実部を示す実部軸、縦軸は積ｈλの値の虚数部を示す虚数部軸である。図１の安定領域曲線１０１、１０２、及び１０３は夫々、前述の（２）式の状態方程式を代表的な数値解法である１次、２次、及び４次のＲｕｎｇｅ−Ｋｕｔｔａ法で解いたときの、発散せずに計算できる安定境界において積ｈλの値をプロットしたものである。積ｈλの計算において、ｈは上記Ｒｕｎｇｅ−Ｋｕｔｔａ法が実行されるときの計算ステップであり、λは（２）式の係数行列Ａと正則な行列Ｐを用いて前述の（３）式によって計算される。

例えば、４次のＲｕｎｇｅ−Ｋｕｔｔａ法は、下記（５）式で示される。

図１で、各次数でのＲｕｎｇｅ−Ｋｕｔｔａ法の計算において、各次数に対応する安定領域曲線１０１、１０２、又は１０３の各内側に積ｈλが存在すれば、対応する次数のＲｕｎｇｅ−Ｋｕｔｔａ法により発散しない計算が可能である。

図１において、破線領域１０４は、Ｓｔｉｆｆな固有値に基づく積ｈλの存在領域である。また、破線領域１０５は、応答を支配する固有値に基づく積ｈλの存在領域である。前述したＳｔｉｆｆ状態方程式を解く場合に、絶対値が大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくして、積ｈλの存在領域１０４が安定領域曲線１０１、１０２、又は１０３の内側に入るようにする必要がある。計算ステップｈを小さくすれば、図１において、積ｈλの存在領域が領域１０４から破線矢印１０６で示される方向に移動して安定領域曲線１０１、１０２、又は１０３の内側に入るため、発散しない計算が可能である。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に基づく積ｈλの存在領域１０５付近では、その計算ステップｈは過剰に小さくなりすぎる。この結果、計算ステップｈを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。

前述した（２）式が前述したＬＴＩシステムの状態方程式となっており、かつ入力が一定値ｕであるとみなせる場合、行列指数関数

を用いて、下記（６）式の一般解を得ることができる。

上記（６）式は、下記（７）式及び（８）式において、Ｎ＝∞（無限大）の場合に相当する。

しかし、Ｎ＝∞（無限大）の場合には、計算が困難である。このため、（７）式及び（８）式において有限のＮを用いると共に、計算誤差を少なくするために微小時間＝ｈの計算ステップで（７）式及び（８）式の差分方程式を計算する手法がある。この場合、Ｘ_ｎはＸ（ｈ×ｎ）の近似値であり、ｕは［ｔ×ｈ，（ｔ＋１）×ｈ］の間一定であるとし、ｕ_ｎとする。上記（７）式及び（８）式の演算法を、本明細書ではＴａｙｌｏｒ展開法と呼ぶ。

図２は、Ｔａｙｌｏｒ展開法において安定に計算できる固有値λと計算ステップｈの積ｈλの複素平面領域を説明する図である。図２において、図１の場合と同様に、横軸は積ｈλの値の実部を示す実部軸、縦軸は積ｈλの値の虚数部を示す虚数部軸である。図２の安定領域曲線２０１、２０２、２０３、２０４、２０５、及び２０６は夫々、前述の（７）式及び（８）式においてN=１、２、４、５、８及び１０とした場合に計算したときの、発散せずに計算できる安定境界において積ｈλの値をプロットしたものである。積ｈλの計算において、ｈは上記（７）式及び（８）式の状態方程式の解法演算が実行されるときの計算ステップｈに等しく、λは（７）式及び（８）式の係数行列Ａと正則な行列Ｐを用いて前述の（３）式によって計算される。

なお、図２において、積ｈλの虚数部の値が−πｉから＋πｉ（ｉは虚数単位）までの間のグレー色に塗りつぶされた範囲２０７は、ナイキスト周波数を超えない範囲であり、その外側の範囲はナイキスト周波数を超える範囲である。積ｈλとナイキスト周波数との関係については、本実施形態の説明において後述する。

また、図２において、図１の場合と同様に、各次数においてでのＴａｙｌｏｒ展開法の計算において、各次数に対応する安定領域曲線２０１、２０２、２０３、２０４、２０５、又は２０６の各内側に積ｈλが存在すれば、発散しない計算が可能である。なお、Ｔａｙｌｏｒ展開法の次数とは前述した（８）式におけるＮの値を指すものである。

なお、（７）式及び（８）式が前述したＬＴＩシステムの状態方程式となっている場合には、前述した（５）式のＲｕｎｇｅ−Ｋｕｔｔａ法は、式を整理すれば４次のＴａｙｌｏｒ展開法となる。４次までのＲｕｎｇｅ−Ｋｕｔｔａ法の安定領域（図１の１０１、１０２、及び１０３）は、Ｔａｙｌｏｒ展開法の安定領域（図２の２０１、２０２、及び２０３）と同一であることが、一般に知られている。

以上説明したＲｕｎｇｅ−Ｋｕｔｔａ法及びＴａｙｌｏｒ展開法は、広く使われる数値計算手法である。しかし、これらの数値計算法は、安定領域に限りがある。そして、Ｓｔｉｆｆ状態方程式の場合、前述したように、不安定になり計算が発散するか、過剰に小さな計算ステップｈを用いなければならないという問題がある。

このような問題を解決するために、複素平面のいわゆる左半平面全体を計算安定領域とする、下記（９）式で示される後退Ｅｕｌｅｒ法又は下記（１０）式で示される台形法（Ｔｕｓｔｉｎ法）と呼ばれる解法が、一般的に用いられている。

しかしながら、（９）式又は（１０）式に示される解法では逆行列を計算する必要がある。一般に逆行列の計算負荷は、行列が“ｎ×ｎ”要素のサイズを有する場合、ｎの３乗に比例すると言われている。従って、ｎが大きなシステムの場合には、計算負荷が極めて高くなる。このため、Ｅｕｌｅｒ法や台形法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。

更に、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次、台形法は２次の近似精度となる。このため、固有値λの絶対値が大きい場合は、これらの解法は、安定動作範囲が広いため発散しないで計算できるものの、計算精度が著しく落ちるというあまり知られていない事実がある。

以上の考察からわかるように、計算装置において絶対値が大きな固有値λを有する方程式を適切に解く演算を実行することは、一般的に困難である。

そこで、本実施形態は、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とすることを目的とする。本実施形態の詳細について、以下に説明する。

本実施形態の動作原理について説明する。まず、一般的な状態方程式である前述の（２）式において、変数ベクトル項Ｘ（ｔ）の係数行列Ａを、対角要素が固有値である対角行列、又はジョルダン標準形行列に変換する。本実施形態において、これらの行列を総称して、固有値包含対角行列と呼ぶことにする。

例えば、（２）式の係数行列Ａ、Ｂが定数であるＬＴＩシステムにおいて、下記（１１）式（前述した（３）式の再掲）に示される変換演算が実行される。

この結果、（２）式の固有値を対角要素に持つ下記（１１）式（前述の（４）式の再掲）に示される対角行列Λが、上述の固有値包含対角行列として算出される。

或いは、（２）式の係数行列Ａのサイズをｎ行ｎ列（ｎは自然数）とし、正則な行列Ｐを用いて、係数行列Ａに対して、下記（１３）式及び（１４）式で示される演算が実行される。この結果、固有値λ_１、・・・、λ_ｓ、・・・、λ_ｍ（１≦ｓ≦ｍ）を対角要素に持つ（１３）式で示されるジョルダン細胞行列と呼ばれる行列を、更に対角要素として持つ（１４）式で示されるジョルダン標準形行列Ｊが算出される。結果的に、ジョルダン標準形行列Ｊは、その対角要素ｓ（１≦ｓ≦ｍ）に固有値λ_ｓを含む、固有値包含対角行列である。

そして、本実施形態では、上記固有値包含対角行列内の各固有値が所定の条件を満たす場合に、その固有値の実部又は虚数部の一方又は両方を所定の定数に変更する。

例えば、対角行列Λ内の固有値λ_ｓ（１≦ｓ≦ｎ）、又はジョルダン標準形行列Ｊ内の固有値λ_ｓ（１≦ｓ≦ｍ）の実部σ_ｓの絶対値｜σ_ｓ｜が所定の実部閾値ａより大きい場合に、上記固有値λ_ｓの実部σ_ｓを所定の実数定数ｂに変更する。より具体的には、任意の数値計算手法、例えば前述の（７）式及び（８）式のＴａｙｌｏｒ展開法の演算に着目する。例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば実数“−３”の絶対値“３”を、上記実部閾値ａとする。そして、固有値λ_ｓの実部σ_ｓの絶対値｜σ_ｓ｜が上記実部閾値ａ＝３より大きい場合に、固有値λ_ｓの実部σ_ｓを所定の実数定数ｂに変更するものである。

また例えば、対角行列Λ内の固有値λ_ｓ（１≦ｓ≦ｎ）、又はジョルダン標準形行列Ｊ内の固有値λ_ｓ（１≦ｓ≦ｍ）の虚数部ω_ｓｉの絶対値｜ω_ｓ｜が所定の虚数部閾値ｃより大きい場合に、上記固有値λ_ｓの虚数部ω_ｓｉを所定の虚数定数ｄｉに変更する。例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば虚数“πｉ”（ｉ：虚数単位）の絶対値“π”を、上記虚数部閾値ｃとする。そして、固有値λ_ｓの虚数部ω_ｓｉの絶対値｜ω_ｓ｜が上記虚数部閾値ｃ＝πより大きい場合に、固有値λ_ｓの虚数部ω_ｓｉを所定の虚数定数ｄｉ（ｄは正の実数）に変更するものである。また例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば虚数“−πｉ”の絶対値“π”を、上記虚数部閾値ｃ′とする。このとき、固有値λ_ｓの虚数部ω_ｓｉの絶対値｜ω_ｓ｜が上記虚数部閾値ｃ′＝πより大きい場合に、固有値λ_ｓの虚数部ω_ｓｉを所定の虚数定数−ｄｉ（ｄは正の実数）に変更するものである。

従って、例えば（１２）式の対角行列Λ内の各対角要素λ_２及びλ_３において、その実部の絶対値｜σ_２｜及び｜σ_３｜が実部閾値ａ＝３より大きく、かつその虚数部の絶対値｜ω_２｜及び｜ω_３｜が虚数部閾値ｃ＝π又はｃ′＝πより大きければ、対角行列Λが下記（１５）式に例示されるように変更されるものである。この行列を、修正対角行列と呼ぶ。修正対角行列は、固有値包含対角行列を変更した修正固有値包含対角行列に対応する。

その後、（２）式に例示される状態方程式が固有値包含対角行列内の固有値を変更して得られる修正固有値包含対角行列に基づく形式に変換される。そして、変換後の状態方程式に基づいて、例えば１０次のＴａｙｌｏｒ展開法により得られる差分方程式に対する解法演算が実行される。

以上の動作原理によれば、（２）式に例示される状態方程式を常に積ｈλに関する安定領域内で解くことが可能となる。このため、状態方程式がたとえＳｔｉｆｆ状態方程式であったとしても、上記安定性が失われることはない。そして、状態方程式が変更された後は、在来のＲｕｎｇｅ−Ｋｕｔｔａ法やＴａｙｌｏｒ展開法により得られる差分方程式に対する解法演算が実行される。この結果、本出願によれば、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、ＨＩＬＳおよびＳＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。

上記動作原理に基づく本発明のいくつかの実施形態の詳細について、以下に説明する。図３は、本発明の第１の実施形態である、ＨＩＬＳ計算装置３００の実施形態を示すブロック図である。ＨＩＬＳ計算装置３００は、物理的な信号入出力手段を有する、いわゆるＨＩＬＳの実施形態である。

ＨＩＬＳ計算装置３００は、特には図示しないが、所謂組込みコンピュータのハードウェアとして、ＣＰＵ（中央演算処理装置）、ＲＯＭ（リードオンリーメモリ）、ＲＡＭ（ランダムアクセスメモリ）、フラッシュメモリ装着装置などを有する。図３に示されるＨＩＬＳ計算装置３００において、通信手段３０１は所謂イーサネットあるいはＣＡＮ通信、またはこれに準ずるものである。固有値計算手段３０２、固有値実部変更手段３０３、固有値虚数部変更手段３０４、行列算出手段３０５、ベクトル算出手段３０８、及び状態方程式計算手段３０９は、上記ＣＰＵが上記フラッシュメモリ装着装置に装着されたメモリカードに記憶されたＨＩＬＳ制御プログラムを上記ＲＡＭに読み出して実行する機能ブロック群である。時計装置３１０は、ＨＩＬＳ計算装置３００のコンピュータが内蔵するタイマ装置のハードウェアである。デジタル信号入力手段３０６、電圧入力手段３０７、電圧出力手段３１１は、上記入出力インタフェース装置及び上記ＲＡＭ又は上記外部記憶装置によって構成される。通信手段３０１には、ローカルエリアネットワークやインターネット等のネットワークまたは専用通信線を介してホストコンピュータ３２０が通信接続される。

デジタル信号入力手段３０６及び電圧入力手段３０７は、計算装置側物理的入力手段として動作し、特には図示しない制御装置の制御装置側物理的出力手段のデジタル信号出力手段及び電圧出力手段に接続される。

デジタル信号入力手段３０６は、制御装置からの複数の信号線の各電圧を監視し、夫々所定の閾値を超えた場合に“１”、超えない場合に“０”の各デジタル入力信号をベクトル算出手段３０８に出力する。

電圧入力手段３０７は、複数のアナログ／デジタル変換器を内蔵し、制御装置からの複数のアナログ入力電圧値をデジタル入力電圧値に変換し、各デジタル入力電圧値をベクトル算出手段３０８に出力する。

ベクトル算出手段３０８は、上記複数のデジタル入力信号の“１”又は“０”値及び複数のデジタル入力電圧値をベクトル要素とする入力ベクトルＵ（前述した（２）式のＵ（ｔ）に対応）を算出する。なお、入力ベクトルＵの第１要素としては値“１”を設定する。この第１要素は、状態方程式が定数項（前述した（２）式のＢＵ（ｔ）＝定数）を必要とした場合に、係数行列Ｂの値により任意の定数を生成するために使用する。

一方、ホストコンピュータ３２０は、一般的なパーソナルコンピュータであってよく、ＨＩＬＳの数学モデルを生成するためのシミュレーションモデリングツールのプログラムを実行することができる。ユーザは、例えば図４に例示されるような機械システムを制御対象とし、その数学モデルを生成したい場合には、上記シミュレーションモデリングツールを実行することにより、例えば図５に例示されるようなブロック線図を、ホストコンピュータ３２０のキーボード、マウス、及びディスプレイを使って作成する。

ホストコンピュータ３２０には更に、例えばシミュレーションモデリングツールの機能として、状態方程式生成ソフトウェアが搭載されている。この状態方程式生成ソフトウェアは、ユーザが図５に例示されるように作成したブロック線図に記載されている積分器の出力値を要素とする、例えば下記（１６）式として例示されるような状態変数ベクトルを生成する。

状態方程式生成ソフトウェアは更に、上記ブロック線図を解析し、積分器に与えられる入力に基づいて、係数行列Ａ及びＢを生成する。例えば、図５に例示されるブロック線図におけるｖ１の出力である第１積分器５０１に接続されている値は、第２積分器５０２及び第３積分器５０３に夫々所定の係数を乗じてｍ１で除算したものである。従って、状態方程式生成ソフトウェアは、下記（１７）として例示される係数行列Ａ及びＢを生成する。

上記のようにして算出される係数行列Ａ及びＢにより、下記（１８）式の状態方程式が構成されることになる。

ここで、Ｆは、制御装置から与えられる入力値である。従って、（１８）式の状態方程式は、前述した（２）式の状態方程式に対応する。

更に、状態方程式生成ソフトウェアは、ユーザが指定することにより、制御装置に接続したい出力値をＹベクトルとして設定することができ、このＹベクトルを生成するための下記（１９）式として示される出力方程式を構成する係数行列Ｃ及びＤを生成する。

なお、図1で示したＳｔｉｆｆな固有値の例（図１の領域１０４を参照）は、図４に例示した機構に対して、以下の諸元を与えたものである。

ｍ１＝３００ ｍ２＝３００ ｋｓ＝１０ ｋ＝１０ ｃ＝１０００

このように、Ｓｔｉｆｆ状態方程式は極めて一般的に存在する。特に、ＳＩＬＳやＨＩＬＳの用途である機械装置や電子回路の模擬においては、いままで避けて通れない難題であった。本実施形態は、これを解決する手段を提供するものである。

以上説明したようにしてホストコンピュータ３２０で生成された係数行列Ａ，Ｂ，Ｃ，及びＤのデータが、ホストコンピュータ３２０からＨＩＬＳ計算装置３００に送られ、図３の通信手段３０１で受信され、記憶される。

図３の説明に戻り、固有値計算手段３０２は、通信手段３０１が受信した上記係数行列Ａを受け取り、この係数行列Ａに基づいて下記の手順により固有値を算出する。

手順１として、下記（２０）式として示される固有方程式を計算する。

係数行列Ａのサイズがｎ行ｎ列の場合、固有値λはｎ次多項式となる。ここで、ベアストウ法として知られた手法を用いることによって、全ての固有値λ_１、・・・、λ_ｓ、・・・、λ_ｎを算出する。これらの固有値λ_ｓ（１≦ｓ≦ｎ）に重複がなければ、これらを対角要素とする対角行列Λを、前述した（１２）式の形で生成する。

次に、固有値λ_ｓ（１≦ｓ≦ｎ）に重複がなければ、手順２として、固有値λ_ｓに対する固有ベクトルｘ_ｓ（１≦ｓ≦ｎ）を、下記（２１）式により算出する。

この計算は、ガウス消去法として知られる連立方程式の解を得るアルゴリズムにより実行される。

固有値λ_ｓ（１≦ｓ≦ｎ）に重複がなければ、更に手順３とし、手順２で算出した固有ベクトルｘ_ｓ（１≦ｓ≦ｎ）を用いて、下記（２２）式として示される正則行列Ｐを構成する。

これと共に、Ｐの逆行列Ｐ^−１を、ガウス・ジョルダン法と呼ばれる手法で計算する。

固有値λ_ｓ（１≦ｓ≦ｎ）に重複がなければ、手順１で算出した対角行列Λを、固有値実部変更手段３０３及び固有値虚数部変更手段３０４に出力する。また、手順３で得られた、正則行列Ｐ及びその逆行列Ｐ^−１と、上記対角行列の逆行列Λ^−１とを、後述するように係数行列Ｂを修正するために、行列算出手段３０５に出力する。

なお、手順１で算出された固有値λ_ｓ（１≦ｓ≦ｎ）に重複がある場合には、固有値計算手段３０２は、対角行列Λの代わりに、前述した（１３）式及び（１４）式で示したジョルダン標準形行列Ｊを算出する。ジョルダン標準形行列Ｊの算出方法は従来知られているが、ここではその説明は省略する。

また、手順１で算出された固有値λ_ｓ（１≦ｓ≦ｎ）に重複がある場合には、重複した固有値に微小値を加減算することにより重複を解除して対角行列Λを作成して、ＳＩＬＳやＨＩＬＳに用いる場合は実用上の問題はない。

次に、固有値実部変更手段３０３及び固有値虚数部変更手段３０４の動作について説明する。図６は、固有値変更手段として一体として動作する固有値実部変更手段３０３及び固有値虚数部変更手段３０４の動作例を示すフローチャートである。

まず、図３の固有値計算手段３０２から、対角行列Λを取得する（図６のステップＳ６０１）。

次に、ｎ個の固有値λ_ｓ（１≦ｓ≦ｎ）について処理を繰り返すための変数ｓに値１をセットする（図６のステップＳ６０２）。

その後、変数ｓの値が上記ステップＳ６０２で値１にセットされた後、図６のステップＳ６０９で変数ｓの値が＋１ずつインクリメントされながら図６のステップＳ６１０で変数ｓの値が値ｎを超えたと判断するまで、図６のステップＳ６０３からステップＳ６０８までの一連の処理を繰り返し実行する。以下、これら一連の処理について順次説明する。

まず、ステップＳ６０１で取得した対角行列Λから、変数ｓの値が示す位置の１つの対角要素の固有値“λ_ｓ＝σ_ｓ＋ω_ｓｉ”を取得する（図６のステップＳ６０３）。ここで、σ_ｓは固有値λ_ｓの実部、ω_ｓｉはｉを虚数単位として固有値λ_ｓの虚数部である。

次に、ステップＳ６０３で取得した固有値λ_ｓの実部σ_ｓの絶対値｜σ_ｓ｜が所定の実部閾値ａより大きいか否かを判定する（図６のステップＳ６０４）。

ステップＳ６０４の判定がＹＥＳならば、固有値λ_ｓの実部σ_ｓを所定の実数定数ｂに変更する（図６のステップＳ６０５）。

ステップＳ６０４の判定がＮＯならば、固有値λ_ｓの実部σ_ｓはそのままとして変更しない。

以上のステップＳ６０４とＳ６０５の処理が、図３の固有値実部変更手段３０３の機能を実現する。

次に、ステップＳ６０３で取得した固有値λ_ｓの虚数部ω_ｓｉの絶対値｜ω_ｓ｜が所定の虚数部閾値ｃより大きいか否かを判定する（図６のステップＳ６０６）。

ステップＳ６０６の判定がＹＥＳならば、固有値λ_ｓの虚数部ω_ｓｉを所定の虚数定数ｄｉ（ｉは虚数単位）に変更する（図６のステップＳ６０７）。

ステップＳ６０６の判定がＮＯならば、固有値λ_ｓの虚数部ω_ｓｉはそのままとして変更しない。

以上のステップＳ６０６とＳ６０７の処理が、図３の固有値虚数部変更手段３０４の機能を実現する。

その後、適宜変更後の“σ_ｓ＋ω_ｓｉ”を、修正対角行列（修正固有値包含対角行列）

において、変数ｓの現在値が示す位置の対角要素の固有値λ_ｓとして設定する（図６のステップＳ６０８）。

ステップＳ６０８の処理の後、変数ｓの値を＋１インクリメントする（図６のステップＳ６０９）。

そして、インクリメント後の変数ｓの値が対角行列Λの次数ｎを超えたか否かを判定する（図６のステップＳ６１０）。

ステップＳ６１０の判定がＮＯならば、ステップＳ６０３の処理に戻って、新たな変数ｓの値が示す対角行列Λ上の対角要素位置の固有値λ_ｓの変更処理を実行する。

ステップＳ６１０の判定がＹＥＳになると、図３の固有値実部変更手段３０３及び固有値虚数部変更手段３０４の今回の計算ステップにおける処理を終了し、修正対角行列

を、後述する図３の行列算出手段３０５及び状態方程式計算手段３０９に出力する。

以上のようにして、本実施形態では、修正対角行列

に基づいて実行される状態方程式が常に積ｈλに関する安定領域内で解かれるように、対角行列Λの各対角要素の固有値を例えば前述した（１５）式に例示される修正対角行列の各対角要素のように変更することが可能となる。

図３の説明に戻り、行列算出手段３０５は、係数行列Ｂに関する以下の修正処理を実行する。まず、前述した（３）式及び（４）式で示されるＬＴＩシステムの状態方程式を前提として、下記（２３）式に示される方程式を設定する。（２３）式の方程式においては、t=∞値（無限大）を、対角行列Λで計算した場合と、修正対角行列

で計算した場合に結果が定常値でありかつ同一であるという条件が設定される。

上記（２３）式より、下記（２４）式が算出される。

行列算出手段３０５は、上記数（２４）式の演算を実行する。即ち、図３において、行列算出手段３０５は、前述した図６のフローチャートに基づいて一体で動作する固有値実部変更手段３０３及び固有値虚数部変更手段３０４から、修正対角行列

を入力する。また、行列算出手段３０５は、固有値計算手段３０２から、前述した正則行列Ｐ、その逆行列Ｐ^−１、及び対角行列の逆行列Λ^−１を入力する。更に、行列算出手段３０５は、通信手段３０１から、係数行列Ｂを入力する。行列算出手段３０５は、これらの入力したデータに基づいて、（２４）式により、係数行列Ｂから修正係数行列

を算出し、状態方程式計算手段３０９に出力する。

最後に、図３の状態方程式計算手段３０９は、前述した（３）式及び（４）式で示されるＬＴＩシステムの状態方程式を前提として、以下に説明する状態方程式の計算処理を実行する。
まず、状態方程式計算手段３０９は、所定の値又はホストコンピュータ３２０からＨＩＬＳ計算装置３００内の通信手段３０１に転送された初期値ベクトル
を用いて、前述した（３）式に基づく下記（２５）式により、初期状態変数ベクトルＺ（０）を算出する。

以後、状態方程式計算手段３０９は、前述の（３）式に基づく下記（２６）式により、状態変数ベクトルＺ（ｔ）を順次算出する。

状態方程式計算手段３０９は、（２６）式の計算を実際には例えば、前述した（７）式及び（８）式を修正して得られる下記（２７）式及び（２８）式による、例えば５次のＴａｙｌｏｒ展開法の差分方程式として実行する。

状態方程式計算手段３０９は、（２７）式及び（２８）式の差分方程式の解法演算を、時計装置３１０が定める所定の時間間隔ｈ秒毎に離散時間ｎをインクリメントしながら実行する。上記（２７）式及び（２８）式の演算において、入力値ｕ_ｎは、ベクトル算出手段３０８から入力する入力ベクトルＵ（ｔ）を、時間間隔ｈでサンプリングして得られるものである。

なお、Ｚ_ｎは一般に複素数であるが、このため、状態方程式計算手段３０９は、（２７）式及び（２８）式により得られる状態変数ベクトルＺ_ｎに対して、下記（２９）式として示される出力方程式を適用することにより、実数である出力ベクトルＹ_ｎを得て電圧出力手段３１１に出力する。

電圧出力手段３１１は、状態方程式計算手段３０９における（２９）式の演算により得られた出力ベクトルＹ_ｎに相当する電圧値を出力する。Ｙ_ｎは時間間隔ｈごとに計算されるため、電圧出力手段３１１の出力は、時間間隔ｈで新しい電圧値に更新される。図３には図示していないが、電圧出力手段３１１内のデジタル／アナログ変換器によりデジタル電圧信号から変換されて出力されたアナログ電圧値は、制御装置の制御装置側物理的入出力手段に入力される。制御装置は、このアナログ電圧値を、制御対象に設置されたセンサ装置からの信号と理解して、制御対象の挙動を判定し、所定の制御出力を出力する。この制御出力は、ＨＩＬＳ計算装置３００のデジタル信号入力手段３０６又は電圧入力手段３０７に入力される。このようにして、制御装置と制御対象であるＨＩＬＳ計算装置３００との間に、フィードバックのクローズドループ（閉ループ）が形成される。

図７は、図３とは異なる本発明の第２の実施形態の構成を示す図であり、ＳＩＬＳ計算装置７００の実施形態を示すブロック図である。図７は、計算装置がＳＩＬＳを構成する例である。図７では、図３の通信手段３０１の代わりに、行列設定手段７０２を有し、状態方程式において用いられる前述したＡ，Ｂ，Ｃ，Ｄ行列が設定される。図７において、図３の場合と同じ参照番号を付したブロックは、図３の場合と同じ機能を有する。

また、ＳＩＬＳ計算装置７００は、制御装置制御則計算手段７０１を備えている。図３のＨＩＬＳ計算装置３００では、制御装置は、デジタル信号入力手段３０６、電圧入力手段３０７、又は電圧出力手段３１１を介してＨＩＬＳ計算装置３００の外部に接続されていた。これに対して、図７では、上記制御装置の動作を、ＳＩＬＳ計算装置７００内の制御装置制御則計算手段７０１によって模擬することができる。

ただし、状態方程式計算手段３０９は計算ステップｈを用いて計算するが、制御装置制御則計算手段７０１の動作は必ずしも物理的にｈ秒毎に実行される必要はない（図７中では「ｈ′」と記載）。この点において、図７のＳＩＬＳ計算装置７００と図３のＨＩＬＳ計算装置３００は異なる。

図７のＳＩＬＳ計算装置７００によれば、それ単体で、内蔵される制御装置制御則計算手段７０１によって制御装置の制御則の動作を検証することが可能である。電圧出力手段３１１の出力は、制御装置制御則計算手段７０１にフィードバックされると共に、ディスプレイ装置７０３に表示させて確認することができる。制御装置制御則計算手段７０１から状態方程式計算手段３０９への入力ベクトルＵ（ｔ）についても、ディスプレイ装置７０３で確認できるようにしてもよい。

図３のＨＩＬＳ計算装置３００又は図７のＳＩＬＳ計算装置７００において、固有値実部変更手段３０３及び固有値虚数部変更手段３０４が固有値変更手段として動作する効果について、以下に説明する。

前述した（２）式が固有値λを持つ場合、その解は前述した（６）式で与えられ、下記（３０）式となる。

λが虚数部を持つ場合は、指数関数は下記（３１）式のオイラー公式より、時間的に振動する波形として現れる。

計算ステップｈ＝１の場合、これはサンプリング周波数＝１秒であるから、ナイキスト周波数は0.5Ｈｚ、周期は２秒であり、これ以上の周波数は表現できない。即ち、図２の領域２０７として示したように、ｈ＝１の計算ステップでは、周波数が0.5Ｈｚ＝πラジアン／秒以上である波形は、次数が高いＴａｙｌｏｒ展開法であれば安定領域に含まれる。この場合、発散しない計算が可能であるが、上記波形は、ｈ＝１のサンプリング周波数では表現されず、ナイキスト周波数で折り返して0.5Ｈｚより低い周波数、所謂エイリアス周波数となって現れる。
従って、固有値λの虚数部の絶対値がπラジアン／秒を超えた場合、この虚数部を例えば±πｉに変更して、その変更後の状態方程式を安定に計算できる５次のＴａｙｌｏｒ展開法で計算してもよい。
この場合、表現できない波形という意味を忠実に再現するならば、０ラジアン／秒としてもよい。
なお、Ｔａｙｌｏｒ展開法としては、前述した（７）式及び（８）式を用いることができる。そして、λ＝πｉの場合、例えば４次の項は、π^４／２４≒４．０５９、５次の項はπ^５／１２０≒２．５５となる。この結果は、Ｎ＝∞において得られる真値に対して、誤差が大きい。
これに対して、例えばＮ＝１０の場合、π^１０／１０！＝０．０２５８となる。この場合、打ち切っても誤差が少ないと判断できる。このため、ｈλ＝πという大きな固有値をある程度の精度で計算するためには、この程度まで次数を上げる必要があり、５次では正しい値は得られない。

図８は、状態方程式が、例えば下記（３２）式の方程式

のように、固有値λ＝±πｉを持つ場合に、異なる計算手法で計算した時系列結果を説明する図である。図８は、ｈ＝１として、台形法、１０次のＴａｙｌｏｒ展開法、及び後退Ｅｕｌｅｒ法の計算結果を、夫々、８０１、８０２、及び８０３として、理論値８０４と共に示している。図８（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図８（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。図８から理解されるように、後退Ｅｕｌｅｒ法の計算結果８０３は、振動波形を表現できない。また、台形法の計算結果８０１は、時間軸対称という特性により２次の割には振動波形に強いが、やはり真値とはかけ離れている。
図８に関する考察より、例えばλ＝±πｉのような固有値は、確かにサンプル時間がｈ＝１の場合は丁度ナイキスト周波数であるため波形を表現し得る。しかし、このような固有値を持つ方程式を計算するためには、計算結果８０４として示される１０次のＴａｙｌｏｒ展開法のような、計算負荷が高い手法を採用せざるを得ない。

そこで、本実施形態では、例えばサンプル時間h=1として、図３又は図７の固有値虚数部変更手段３０４において、前述した虚数部閾値ｃがπ／２（＝０．２５Ｈｚ）に設定される。そして、固有値の虚数部の絶対値が上記虚数部閾値ｃ＝π／２より大きなものについては、その固有値の虚数部が全て所定の虚数定数ｄｉ＝π／２ｉ又はｄｉ＝−π／２ｉに変更される。
図９は、固有値λ＝±π／２ｉを異なる計算手法で計算した時系列結果を説明する図である。図８の場合と同様に、図９（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図９（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。Ｔａｙｌｏｒ展開法の４次の項は、（π／２）^４・２４＝０．２５４となる。このため、４次で打ち切られるＴａｙｌｏｒ展開法、又は４次のＲｕｎｇｅ−Ｋｕｔｔａ法を用いても、図９に示されるように、十分な計算精度を得ることができる。なお、図１及び図２に関して前述したように、前述した（２）式についてＬＴＩが仮定でき（３）式及び（４）式のように表現できるならば、４次以下の次数については、Ｔａｙｌｏｒ展開法とＲｕｎｇｅ−Ｋｕｔｔａ法は同一である。
この場合、図９のＲｕｎｇｅ−Ｋｕｔｔａ法の計算結果９０１は、本来の方程式が与える振動波形の周波数の理論値８０４とは明らかに異なる。しかし、従来用いられている代表的な手法である台形法の計算結果８０１であっても、本来の周波数の波形を計算することはできない。このため、少なくとも振動現象を表現できるという点では大差はなく、固有値λを±π／２ｉに変換して、Ｒｕｎｇｅ−Ｋｕｔｔａ法により計算した手法は計算効率に勝る。

以上のように、そもそもにおいて例えばサンプル時間h=1の場合は固有値の虚数部の絶対値がπより大きければ、この固有値による応答は表現することができない。また、固有値の虚数部の絶対値がπ以下の値であっても、１０次のＴａｙｌｏｒ展開法のような特殊な手法を除き、実用になる数値計算法が存在せず、正確な波形再現はできない。そこで、本実施形態では、正確な計算が可能な位置まで固有値を変更することで、結果的により適切な計算結果を得ることがものである。

次に固有値の実部の影響について説明する。図１０は、指数関数の時間応答を異なる計算手法で計算した結果を説明する図である。横軸は時間（秒）、縦軸はｘの計算結果である。図１０は、指数関数が負の実部を持つ場合であって、ｔ＝０での初期値を１とした場合の時間応答である。
ｅ^λｔの最初の１ステップ、即ちｔ＝ｈの値は、

λ＝−１ ｈ＝１：ｅ^−ｔ ＝０．３６８
λ＝−２ ｈ＝１：ｅ^−２ｔ＝０．１３５
λ＝−３ ｈ＝１：ｅ^−３ｔ＝０．０５０
λ＝−４ ｈ＝１：ｅ^−４ｔ＝０．０１８３

である。従って、例えば固有値がλ＝−３より小さければ、その固有値の時間応答特性がλ＝−３のときの時間応答特性と同じである、として実用上は問題はない。２ステップ以降（ｔ＝２ｈ，３ｈ，４ｈ）では、時間応答特性は更に０に近似する。従って、固有値がλ＝−３より小さい場合と固有値がλ＝−３の場合とで、それらの差異は全くなくなる。
なお、“λ＝−３”のような原点から離れた固有値は、図６に示すように台形積分法でも大きな誤差がのこり、正確に計算するなら、例えばＴａｙｌｏｒ展開法の１０次のような計算アルゴリズム（図７では分かりやすいように８次を使用）を必要とするなど、虚数部において説明したように、そもそもにおいて固有値の絶対値が大きい場合は、前述の（７）式及び（８）式において次数Ｎを相当に大きくしないと収束しないため、台形法のような次数の小さな計算アルゴリズムでは、打切り誤差が大きく、正確な計算ができない。
そのため、“λ＝−３”より小さい固有値を“λ＝−３”に置き換えたことによる差異は、“λ＝−３”のような絶対値の大きな固有値を持つ状態方程式を、台形法のような現在広く知られている計算手法で計算した場合に避け得ない計算誤差を考慮すれば、無視できるものである。更に、“−２”程度を閾値にして、これより小さな固有値をすべて“λ＝−２”に置き換えて、４次のＲｕｎｇｅ−Ｋｕｔｔａ法を使用して“λ＝−２”のときの時間応答を正確に計算したほうが、結果的に妥当な計算結果を得ることができる場合もある。

なお広く知られたことであるが、固有値が、下記（３３）式で示される任意の複素数である場合を考える。

この場合、下記（３４）式が成立する。

従って、複素数の固有値に対する時間応答特性は、図８（ｂ）又は図９（ｂ）に例示した虚数部の時間応答特性と、図１０に例示した実部の時間応答特性との積となり、図１１に例示されるような特性になる。従って、固有値の実部の時間応答特性と虚数部の時間応答特性は、独立に考えることができる。

次に、本発明の第３の実施形態について説明する。この実施形態は、図３の本発明の第１の実施形態と同様に、ＨＩＬＳ計算装置の実施形態である。図１２は、本発明の第３の実施形態であるＨＩＬＳ計算装置１２００の実施形態を示すブロック図である。図１２において、図３の場合と同じ参照番号を付したブロックは、図３の場合と同じ機能を有する。図１２のＨＩＬＳ計算装置１２００の構成が図３のＨＩＬＳ計算装置３００の構成と異なる点は、まず、図３の固有値実部変更手段３０３及び固有値虚数部変更手段３０４の代わりに、図１２の固有値変更手段１２０１が具備される点である。次に、図１２の行列算出手段１２０２の動作が、図３の行列算出手段３０５の動作とは異なる。更に、図１２の状態方程式計算手段１２０３の動作が図３の状態方程式計算手段３０９の動作と異なる。また、補足的に、図１２のＨＩＬＳ計算装置１２００には、制御装置１２１０が接続される。制御装置１２１０内の特には図示しない第１の制御装置側物理的出力手段及び第２の制御装置側物理的出力手段は、夫々、ＨＩＬＳ計算装置１２００内のデジタル信号入力手段３０６及び電圧入力手段３０７に接続される。デジタル信号入力手段３０６及び電圧入力手段３０７は、計算装置側物理的入力手段として機能する。一方、ＨＩＬＳ計算装置１２００内の電圧出力手段３１１は、制御装置１２１０内の特には図示しない制御装置側物理的入力手段に接続される。電圧出力手段３１１は、計算装置側物理的出力手段として機能する。

まず、図１２の固有値変更手段１２０１の詳細動作について説明する。前述した図３の本発明の第１の実施形態や図７の本発明の第２の実施形態においては、固有値λの実部σと虚数部ωｉは、夫々、固有値実部変更手段３０３及び固有値虚数部変更手段３０４によって独立に変更可能であった。ここで、Ｒｕｎｇｅ−Ｋｕｔｔａ法やＴａｙｌｏｒ展開法では、図１や図２に示されるように、積λｈの安定領域が概ね半円を描く。従って、これらの数値計算法を使用する場合には、使用する計算アルゴリズムによって、固有値λの実部σ及び虚数部ωｉの間になんらかの拘束条件を設けることが、適切と考えられる。そこで、固有値変更手段１２０１では、使用する数値計算アルゴリズムによって定められる所定の値γ及びμを用いて、実部σ及び虚数部ωの線形結合、即ち下記（３６）式及び（３８）式を満たす関係が設定される。

即ち、固有値変更手段１２０１は、以下の動作を実行する。まず、虚数部ωｉの係数ωが０以上の場合、実部σ及び虚数部ωｉが、所定の値γ及びμを用いて、下記（３５）式の関係を満たすか否かが判定される。

上記（３５）の関係が満たされるならば、下記（３６）式の関係が成立する。

上記（３５）式及び（３６）式の関係より、図１２の固有値変更手段１２０１は、固有値λの虚数部ωｉの係数ωが０以上である場合において、（３５）式の関係が満たされる場合には、（３６）式の関係が成立するようにσ又はωの一方又は両方を変更して固有値λを変更する。
一方、虚数部ωｉの係数ωが０未満の場合、実部σ及び虚数部ωｉが、所定の値γ及びμを用いて、下記（３７）式の関係を満たすか否かが判定される。

上記（３７）の関係が満たされるならば、下記（３８）式の関係が成立する。

上記（３７）式及び（３８）式の関係より、図１２の固有値変更手段１２０１は、固有値λの虚数部ωｉの係数ωが０未満である場合において、（３７）式の関係が満たされる場合には、（３８）式の関係が成立するようにσ又はωの一方又は両方を変更して固有値λを変更する。

図１３は、第３の実施形態における固有値λの変更方法を説明する図である。図１３に示される複素平面図は、第１の実施形態における図２の図と同じである。図１３の安定領域曲線１３０１、１３０２、１３０３、１３０４、１３０５、及び１３０６は夫々、図２の安定領域曲線２０１、２０２、２０３、２０４、２０５、及び２０６と同じである。

今、前述した（３６）式は、図１３の複素平面上では、直線１３１０に対応する。また、前述した（３５）式が満たされる領域は、例えば、虚数部の係数ωが図１３の縦軸上で０以上の領域内の上記安定領域曲線１３０１〜１３０６の外側の領域である。例えば座標１３０２が、上記領域に相当する。従って、σとωが、座標１３０２で代表される関係にあるときには、σとωの座標が例えば直線１３１０上の位置１３２１の座標に移動するように、σ又はωの一方の値或いは両方の値が変更される。これにより、変更後の固有値は、安定領域に入ることになる。

同様に、前述した（３８）式は、図１３の複素平面上では、直線１３１１に対応する。また、前述した（３７）式が満たされる領域は、例えば、虚数部の係数ωが図１３の縦軸上で０未満の領域内の上記安定領域曲線１３０１〜１３０６の外側の領域である。例えば座標１３２２が、上記領域に相当する。従って、σとωが、座標１３２２で代表される関係にあるときには、σとωの座標が例えば直線１３１１上の位置１３２３の座標に移動するように、σ又はωの一方の値或いは両方の値が変更される。これにより、変更後の固有値は、安定領域に入ることになる。

上記固有値の変更原理に基づく、固有値変更手段１２０１の動作について説明する。図１４は、固有値変更手段１２０１の動作例を示すフローチャートである。

まず、図３の固有値計算手段３０２から、対角行列Λを取得する（図１４のステップＳ１４０１）。

次に、ｎ個の固有値λ_ｓ（１≦ｓ≦ｎ）について処理を繰り返すための変数ｓに値１をセットする（図１４のステップＳ１４０２）。

その後、変数ｓの値が上記ステップＳ１４０２で値１にセットされた後、図１４のステップＳ１４１０で変数ｓの値が＋１ずつインクリメントされながら図１４のステップＳ１４１１で変数ｓの値が値ｎを超えたと判断するまで、図１４のステップＳ１４０３からステップＳ１４０９までの一連の処理を繰り返し実行する。以下、これら一連の処理について順次説明する。

まず、ステップＳ１４０１で取得した対角行列Λから、変数ｓの値が示す位置の１つの対角要素の固有値“λ_ｓ＝σ_ｓ＋ω_ｓｉ”を取得する（図１４のステップＳ１４０３）。ここで、σ_ｓは固有値λ_ｓの実部、ω_ｓｉはｉを虚数単位として固有値λ_ｓの虚数部である。

次に、ステップＳ１４０３で取得した固有値λ_ｓの虚数部ω_ｓの値が０以上であるか否かを判定する（図１４のステップＳ１４０４）。

ステップＳ１４０４の判定がＹＥＳならば、前述した（３５）式に対応する下記（３９）式が満たされるか否かを判定する（図１４のステップＳ１４０５）。

ステップＳ１４０５の判定がＹＥＳならば、前述した（３６）式に対応する下記（４０）式の関係が成立するように、例えば下記（４１）式に従って実部σ_ｓの値を変更する（図１４のステップＳ１４０６）。

なお、（４０）式の関係が成立するように、実部σ_ｓと虚数部ω_ｓの両方が変更されてもよい。或いは、（４０）式の関係が成立するように、虚数部ω_ｓのみが変更されてもよい。

ステップＳ１４０５の判定がＮＯならば、上記ステップＳ１４０６の変更処理は実行されず、実部σ_ｓと虚数部ω_ｓは元の値のままとされる。

前述したステップＳ１４０４の判定がＮＯならば、前述した（３７）式に対応する下記（４２）式が満たされるか否かを判定する（図１４のステップＳ１４０７）。

ステップＳ１４０７の判定がＹＥＳならば、前述した（３８）式に対応する下記（４３）式の関係が成立するように、例えば下記（４４）式に従って実部σ_ｓの値を変更する（図１４のステップＳ１４０８）。

なお、（４３）式の関係が成立するように、実部σ_ｓと虚数部ω_ｓの両方が変更されてもよい。或いは、（４３）式の関係が成立するように、虚数部ω_ｓのみが変更されてもよい。

ステップＳ１４０７の判定がＮＯならば、上記ステップＳ１４０８の変更処理は実行されず、実部σ_ｓと虚数部ω_ｓは元の値のままとされる。

以上のステップＳ１４０６の処理の後、ステップＳ１４０５の判定がＮＯの場合、ステップＳ１４０８の処理の後、又はステップＳ１４０７の判定がＮＯの場合に、適宜変更後の“σ_ｓ＋ω_ｓｉ”を、修正対角行列（修正固有値包含対角行列）

において、変数ｓの現在値が示す位置の対角要素の固有値λ_ｓとして設定する（図１４のステップＳ１４０９）。

ステップＳ６０８の処理の後、変数ｓの値を＋１インクリメントする（図１４のステップＳ１４１０）。

そして、インクリメント後の変数ｓの値が対角行列Λの次数ｎを超えたか否かを判定する（図１４のステップＳ１４１１）。

ステップＳ１４１１の判定がＮＯならば、ステップＳ１４０３の処理に戻って、新たな変数ｓの値が示す対角行列Λの対角要素位置の固有値λ_ｓの変更処理を実行する。

ステップＳ１４１１の判定がＹＥＳになると、図１２の固有値変更手段１２０１の今回の計算ステップにおける処理を終了し、修正対角行列

を、後述する図１２の行列算出手段３０５に出力する。

以上のようにして、本実施形態では、実行される状態方程式が常に積ｈλに関する安定領域内で解かれるように、対角行列Λの各対角要素の固有値を変更することが可能となる。

図１２の説明に戻り、次に行列算出手段１２０２について説明する。図３の第１の実施形態では、状態方程式計算手段３０９は、複素固有値を対角要素に持つ対角行列Λを用いて、前述した（３）式及び（４）式に基づいて、解法演算を実行した。しかし、この解法演算は複素演算となるため、対角行列Λの特性によっては、計算時間を要する場合があり得る。

そこで、第３の実施形態では、行列算出手段１２０２が、固有値変更手段１２０１で算出された修正対角行列

から、修正係数行列

及び

を得る。そして、図１２の状態方程式計算手段１２０３が、これらの修正係数行列に基づいて修正した状態方程式に対して、実数演算による計算を実行できるようにする。

まず、対角行列Λを得るために使用した前述した（３）式における正則行列Ｐ及びその逆行列Ｐ^−１に基づいて、下記（４５）式の変換式が得られる。

行列算出手段１２０２は、上記（４５）式の演算を実行する。即ち、図１２において、行列算出手段１２０２は、前述した図１４のフローチャートに基づいて動作する固有値変更手段１２０１から、修正対角行列

を入力する。また、行列算出手段１２０２は、固有値計算手段３０２から、前述した正則行列Ｐ及びその逆行列Ｐ^−１を入力する。行列算出手段１２０２は、これらの入力したデータに基づいて、（４５）式により、係数行列Ａに対応する修正係数行列

を算出し、状態方程式計算手段１２０３に出力する。
なお、行列算出手段１２０２は、状態方程式計算手段１２０３における後述する状態方程式の計算のための補足として、上記修正係数行列の逆行列

も算出し、状態方程式計算手段１２０３に出力する。

また、行列算出手段１２０２は、係数行列Ｂについても修正処理を実行する。図３の第１の実施形態において行列算出手段３０５は、前述した（２４）式の演算により、係数行列Ｂを修正係数行列

に修正した。この（２３）式と、前述した（３）式、及び上記（４５）式に基づいて、下記（４６）式の演算式が得られる。

行列算出手段１２０２は、上記（４６）式の演算を実行する。即ち、図１２において、行列算出手段１２０２は、通信手段３０１から係数行列Ａを入力し、それに対応する逆行列Ａ^−１を算出する。また、行列算出手段１２０２は、通信手段３０１から係数行列Ｂを入力する。行列算出手段１２０２は、これらのデータと、前述の（４５）式により算出した修正係数行列

とに基づいて、（４６）式により、修正係数行列

を算出し、状態方程式計算手段３０９に出力する。

最後に、図３の状態方程式計算手段３０９は、下記（４７）式に示される状態方程式を設定する。

状態方程式計算手段３０９は、（４７）式の計算を実際には例えば、前述した（７）式及び（８）式を修正して得られる下記（４８）式及び（４９）式による、例えば５次のＴａｙｌｏｒ展開法の差分方程式として実行する。

このようにして、図１２から図１４で説明した第３の実施形態では、固有値変更手段１２０１での固有値の変更を更に適切に実行できると共に、状態方程式計算手段１２０３での状態方程式の解法演算を更に効率的に実行することが可能となる。

次に、第４の実施形態について説明する。前述した（２）式の計算方法として、例えば前述した（２７）式及び（２８）式においてＮ＝４として４次のＴａｙｌｏｒ展開法を用いた場合を仮定する。この場合、固有値包含対角行列または修正固有値包含対角行列が含む或る固有値λsと或る計算ステップｈに対しては、下記（５０）式及び（５１）式として示される計算が実質的に計算される。

この時α＝０となる積ｈλ_ｓが存在する。具体的には、ｉを虚数単位として、下記（５２）式又は（５３）式の場合に、上記（５１）式においてα＝０となる。

の場合に、上記（５１）式において、α＝０となる。この条件を満たすhλｓをゼロ点と言うが、これらのゼロ点を複素平面上で表すと、図１５に示されるごとくとなる。図１５に示される複素平面図は、第１の実施形態における図２の図と同じである。図１５の安定領域曲線１５０１、１５０２、１５０３、及び１５０４は夫々、図２の安定領域曲線２０１、２０２、２０３、及び２０４と同じである。図１５において、座標点１５１０及び１５１１が（５２）式の場合に対応し、座標点１５１２及び１５１３が（５３）式の場合に対応する。図１５から理解されるように、これらのゼロ点は、４次のＴａｙｌｏｒ展開法の安定領域曲線１５００３の内側、即ち安定領域内に存在することがわかる。

このような固有値の場合、前述した（５０）式は、下記（５４）式のようになる。

この値はｕ_ｎが一定であれば、ｘ_∞の値、つまり最終値であり、最初の１ステップで最終値に収束してしまい、過渡的な値を得ることができない。最終結果を得るだけの計算であればよいが、ＳＩＬＳやＨＩＬＳの場合過渡応答も極めて重要であるため、この状態は好ましいものではない。

そこで、第４の実施形態では、このような固有値は、例え計算アルゴリズムの安定領域内であっても、他の値に変更する。第４の実施形態では、前述した図１２の第３の実施形態の場合と同様の構成において、固有値変更手段１２０１が、（５１）式においてα＝０となる条件を充足することを判別するゼロ点検出手段を更に備える。そして、このゼロ点検出手段は、上記条件が充足されたことを判別した場合に、固有値λ_ｓを他の固有値に変更する、例えば、下記（５５）式のように固有値の変更を行う。あるいは、５次のＴａｙｌｏｒ展開法で計算する。

前述した（２）式において、下記（５６）式の場合には、その固有値は前述した（５２）式の値となり、前述した（５１）式に対する判別はα＝０となる。

図１６は、係数行列Ａ、Ｂを持つ前述した（２）式を、下記の３つの条件で計算した場合における、状態ベクトルＸの第１要素ｘ１と第２要素ｘ２の時系列計算結果を示す図である。

１６０１：オリジナルの４次のＴａｙｌｏｒ展開法で計算した場合

１６０２：前述の（５５）式のように固有値の変更を行った上で、
４次のＴａｙｌｏｒ展開法で計算した場合

１６０３：５次のＴａｙｌｏｒ展開法で計算した場合

また、図１６には、理論値のオリジナルの時系列計算結果１６０４と、理論値×０．７５の時系列計算結果１６０５も、一緒に示されている。
図１６では、２次方程式

の固有値が、図１２の第３の実施形態の場合と同様の固有値変更手段１２０１により、前述の（５２）式の場合の値とされる。
次に、図１２の第３の実施形態の場合と同様の行列算出手段１２０２により、修正係数行列

及び

が生成される。
そして、図１２の第３の実施形態の場合と同様の状態方程式計算手段１２０３で前述の（４７）式の状態方程式が計算される。図１６（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図１６（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。

図１６において、固有値が前述した（５２）式として示される値のままの計算結果１６０１では、振動現象は表現できていない。これに対して、第４の実施形態において固有値が０．７５倍された計算結果１６０２では、真値１６０４とは異なるものの、振動現象は良く表現されていることがわかる。５次のＴａｙｌｏｒ展開法の計算結果１６０３も、やはり真値１６０４とは異なる。
計算結果が真値と異なるのは、固有値の絶対値が大きく４次や５次では計算誤差が大きいためであるので、例えば１０次のＴａｙｌｏｒ展開法などを使用すれば、真値に近い値を得ることができる。
固有値が前述した（５３）式として示される値も、第４の実施形態は同様に適用できる。

上述の第４の実施形態では、固有値変更手段１２０１内のゼロ点検出手段が条件が充足されたことを判別した場合に、例えば下記（５５）式のように、固有値λ_ｓを他の固有値に変更した。これに対し、ゼロ点検出手段が条件が充足されたことを判別した場合に、固有値を変更するのではなく、計算次数Ｎを変更するようにしてもよい。

以上説明したように、固有値計算手段３０２と、固有値実部変更手段３０３／固有値虚数部変更手段３０４又は固有値変更手段１２０１とを有する上述の各実施形態によるＨＩＬＳ，ＳＩＬＳ計算装置は、いわゆるＳｔｉｆｆ状態方程式であっても、計算精度が高く、計算速度の速い数値計算手法を用いることが可能となる。
これらの構成がＳＩＬＳに適用された場合は、計算精度が向上すると同時に計算時間を短縮することが可能となる。
上記構成がＨＩＬＳに適用された場合には、従来では不可能であったＳｔｉｆｆ状態方程式で表される制御対象を、リアルタイム処理として模擬することが可能となる。
この結果、制御装置の開発を飛躍的に効率化することが可能となり、産業界の発展に寄与するものである。

１０１、１０２、１０３、２０１、２０２、２０３，２０４、２０５、２０６、１３０１、１３０２、１３０３、１３０４、１３０５、１３０６、１５０１、１５０２、１５０３、１５０４、 安定領域曲線
２０７ ナイキスト周波数を超えない範囲
３００、１２００ ＨＩＬＳ計算装置
３０１ 通信手段
３０２ 固有値計算手段
３０３ 固有値実部変更手段
３０４ 固有値虚数部変更手段
３０５ 行列算出手段
３０６ デジタル信号入力手段
３０７ 電圧入力手段
３０８ ベクトル算出手段
３０９ 状態方程式計算手段
３１０ 時計装置
３１１ 電圧出力手段
３２０ ホストコンピュータ
７００ ＳＩＬＳ計算装置
７０１ 制御装置制御則計算手段
７０２ 行列設定手段
７０３ ディスプレイ装置

Claims (13)

1. 制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置であって、
   制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間において不変な定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定手段と、
   前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算手段と、
   前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部の一方又は両方を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更手段と、
   前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換手段と、
   前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算手段と、
   を備える制御対象の動作を模擬する計算装置。
2. 前記計算装置はハードウェア・イン・ザ・ループ・シミュレータ装置として動作し、
   前記計算装置は、前記制御装置が備える制御装置側物理的出力手段に接続される計算装置側物理的入力手段と、前記制御装置が備える制御装置側物理的入力手段に接続される計算装置側物理的出力手段とを更に備え、
   前記入力ベクトル項の各要素として、前記制御装置側物理的出力手段から前記計算装置側物理的入力手段に入力する電圧値又は電流値が設定され、
   前記変数ベクトル項の各要素として、前記計算装置側物理的出力手段から前記制御装置側物理的入力手段へ出力される電圧値または電流値、あるいは電磁光学的入力手段の検知対象となる物理値、あるいは通信入力装置に対する通信出力が計算される、
   請求項１に記載の制御対象の動作を模擬する計算装置。
3. 前記固有値変更手段は、前記固有値包含対角行列内の各前記固有値λ_ｓの実部σ_ｓの絶対値｜σ_ｓ｜が所定の実部閾値ａより大きい場合に、前記固有値λ_ｓの実部σ_ｓを所定の実数定数ｂに変更する、請求項１又は２に記載の制御対象の動作を模擬する計算装置。
4. 前記固有値変更手段は、ｉを虚数単位として、前記固有値包含対角行列内の各前記固有値λ_ｓの虚数部ω_ｓｉの絶対値｜ω_ｓ｜が所定の虚数部閾値ｃより大きい場合に、前記固有値λ_ｓの虚数部ω_ｓを所定の虚数定数ｄｉに変更する、請求項１又は２に記載の制御対象の動作を模擬する計算装置。
5. 前記固有値変更手段は、ｉを虚数単位として、前記固有値包含対角行列内の任意の固有値λ_ｓが、
   

   

   で示される複素固有値であるとき、μ及びγを夫々所定の定数として、
   

   

   で示される関係がある場合に、前記（２）式の関係を満たすようにσ_ｓ又はω_ｓｉの一方又は両方を変更して前記固有値λ_ｓを変更する、
   請求項１又は２に記載の制御対象の動作を模擬する計算装置。
6. 前記方程式設定手段は、
   

   

   を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
   

   

   を設定し、
   前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
   

   

   で示される演算により、固有値λ_１、・・・、λ_ｓ、・・・、λ_ｎ（１≦ｓ≦ｎ）を前記対角要素とする前記固有値包含対角行列である対角行列Λに変換し、
   前記固有値変更手段は、前記（４）式で示される前記対角行列Λ内の各前記固有値λ_ｓ＝σ_ｓ＋ω_ｓｉ（１≦ｓ≦ｎ）が所定の条件を満たす場合に、前記固有値λ_ｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
   

   

   を出力し、
   前記方程式変換手段は、前記係数行列Ｂを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記（３）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
   

   

   

   として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
   

   

   で示される出力方程式を出力し、
   前記状態方程式計算手段は、前記（６）式と（７）式とで示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｚ（ｔ）を計算し、更に、前記（８）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
   請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
7. 前記方程式設定手段は、
   

   

   を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
   

   

   を設定し、
   前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
   

   

   で示される演算により、固有値λ_１、・・・、λ_ｓ、・・・、λ_ｎ（１≦ｓ≦ｎ）を前記対角要素とする前記固有値包含対角行列である対角行列Λに変換し、
   前記固有値変更手段は、前記（１０）式で示される前記対角行列Λ内の各前記固有値λ_ｓ＝σ_ｓ＋ω_ｓｉ（１≦ｓ≦ｎ）が所定の条件を満たす場合に、前記固有値λ_ｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
   

   

   を出力し、
   前記方程式変換手段は、前記係数行列Ａを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記係数行列Ｂを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記（９）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
   

   

   として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
   

   

   で示される出力方程式を出力し、
   前記状態方程式計算手段は、前記（１３）式で示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｘ（ｔ）を計算し、更に、前記（１４）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
   請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
8. 前記方程式設定手段は、
   

   

   を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
   

   

   を設定し、
   前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
   

   

   

   で示される演算により、固有値λ_１、・・・、λ_ｓ、・・・、λ_ｍ（１≦ｓ≦ｍ）を前記対角要素に含む前記固有値包含対角行列であるジョルダン標準形行列Ｊに変換し、
   前記固有値変更手段は、前記（１６）式及び前記（１７）式で示される前記ジョルダン標準形行列Ｊ内の各前記固有値λ_ｓ＝σ_ｓ＋ω_ｓｉ（１≦ｓ≦ｍ）が所定の条件を満たす場合に、前記固有値λ_ｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
   

   

   を出力し、
   前記方程式変換手段は、前記係数行列Ｂを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記（１５）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
   

   

   

   として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
   

   

   で示される出力方程式を出力し、
   前記状態方程式計算手段は、前記（１９）式と（２０）式とで示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｚ（ｔ）を計算し、更に、前記（２１）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
   請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
9. 前記方程式設定手段は、
   

   

   を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
   

   

   を設定し、
   前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
   

   

   

   で示される演算により、固有値λ_１、・・・、λ_ｓ、・・・、λ_ｍ（１≦ｓ≦ｍ）を前記対角要素に含む前記固有値包含対角行列であるジョルダン標準形行列Ｊに変換し、
   前記固有値変更手段は、前記（２３）式及び前記（２４）式で示される前記ジョルダン標準形行列Ｊ内の各前記固有値λ_ｓ＝σ_ｓ＋ω_ｓｉ（１≦ｓ≦ｍ）が所定の条件を満たす場合に、前記固有値λ_ｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
   

   

   を出力し、
   前記方程式変換手段は、前記係数行列Ａを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記係数行列Ｂを、
   

   

   で示される演算により、修正係数行列
   

   

   に変換し、前記（２２）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
   

   

   として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
   

   

   で示される出力方程式を出力し、
   前記状態方程式計算手段は、前記（２７）式で示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｘ（ｔ）を計算し、更に、前記（２８）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
   請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
10. 前記状態方程式計算手段は、Ｎを任意の自然数として、
    

    

    

    で示される前記差分方程式に対する、Ｔａｙｌｏｒ展開法又はこれと等価な計算式となるＲｕｎｇｅ−Ｋｕｔｔａ法による解法演算を実行し、
    或る計算ステップｈに対して、前記固有値λ_ｓが、
    

    

    で示される条件を充足することを判別するゼロ点検出手段を更に備え、
    前記固有値変更手段は、前記ゼロ点検出手段が前記（３１）式で示される条件が充足されたことを判別した場合に、前記固有値λ_ｓを他の固有値に変更する、
    請求項１、２、３、４、５、７、又は９の何れか一項に記載の制御対象の動作を模擬する計算装置。
11. 前記状態方程式計算手段は、Ｎを任意の自然数として、
    

    

    

    で示される前記差分方程式に対する、Ｔａｙｌｏｒ展開法又はこれと等価な計算式となるＲｕｎｇｅ−Ｋｕｔｔａ法による解法演算を実行し、
    或る計算ステップｈに対して、前記固有値λ_ｓが、
    

    

    で示される条件を充足することを判別するゼロ点検出手段を更に備え、
    前記方程式変換手段は、前記ゼロ点検出手段が前記（３５）式で示される条件が充足されたことを判別した場合に、前記Ｎの値を変更する、
    請求項１、２、３、４、５、７、又は９の何れか一項に記載の制御対象の動作を模擬する計算装置。
12. 制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置に、
    制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定処理と、
    前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算処理と、
    前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更処理と、
    前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換処理と、
    前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算処理と、
    を実行させる、計算装置において状態方程式を解く演算の実行方法。
13. 制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置に、
    制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定処理と、
    前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算処理と、
    前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更処理と、
    前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換処理と、
    前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算処理と、
    を実行させるためのプログラム。

Images