(57)【特許請求の範囲】
【請求項１】
制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置であって、
制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間において不変な定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定手段と、
前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算手段と、
前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部の一方又は両方を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更手段と、
前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換手段と、
前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算手段と、
を備える制御対象の動作を模擬する計算装置。
【請求項２】
前記計算装置はハードウェア・イン・ザ・ループ・シミュレータ装置として動作し、
前記計算装置は、前記制御装置が備える制御装置側物理的出力手段に接続される計算装置側物理的入力手段と、前記制御装置が備える制御装置側物理的入力手段に接続される計算装置側物理的出力手段とを更に備え、
前記入力ベクトル項の各要素として、前記制御装置側物理的出力手段から前記計算装置側物理的入力手段に入力する電圧値又は電流値が設定され、
前記変数ベクトル項の各要素として、前記計算装置側物理的出力手段から前記制御装置側物理的入力手段へ出力される電圧値または電流値、あるいは電磁光学的入力手段の検知対象となる物理値、あるいは通信入力装置に対する通信出力が計算される、
請求項１に記載の制御対象の動作を模擬する計算装置。
【請求項３】
前記固有値変更手段は、前記固有値包含対角行列内の各前記固有値λｓの実部σｓの絶対値｜σｓ｜が所定の実部閾値ａより大きい場合に、前記固有値λｓの実部σｓを所定の実数定数ｂに変更する、請求項１又は２に記載の制御対象の動作を模擬する計算装置。
【請求項４】
前記固有値変更手段は、ｉを虚数単位として、前記固有値包含対角行列内の各前記固有値λｓの虚数部ωｓｉの絶対値｜ωｓ｜が所定の虚数部閾値ｃより大きい場合に、前記固有値λｓの虚数部ωｓを所定の虚数定数ｄｉに変更する、請求項１又は２に記載の制御対象の動作を模擬する計算装置。
【請求項５】
前記固有値変更手段は、ｉを虚数単位として、前記固有値包含対角行列内の任意の固有値λｓが、
【数１】
0006905296.tif 000080
で示される複素固有値であるとき、μ及びγを夫々所定の定数として、
【数２】
0006905296.tif 000081
で示される関係がある場合に、前記（２）式の関係を満たすようにσｓ又はωｓｉの一方又は両方を変更して前記固有値λｓを変更する、
請求項１又は２に記載の制御対象の動作を模擬する計算装置。
【請求項６】
前記方程式設定手段は、
0006905296.tif 000082
を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
【数３】
0006905296.tif 000083
を設定し、
前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
【数４】
0006905296.tif 000084
で示される演算により、固有値λ１、・・・、λｓ、・・・、λｎ（１≦ｓ≦ｎ）を前記対角要素とする前記固有値包含対角行列である対角行列Λに変換し、
前記固有値変更手段は、前記（４）式で示される前記対角行列Λ内の各前記固有値λｓ＝σｓ＋ωｓｉ（１≦ｓ≦ｎ）が所定の条件を満たす場合に、前記固有値λｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
0006905296.tif 000085
を出力し、
前記方程式変換手段は、前記係数行列Ｂを、
【数５】
0006905296.tif 000086
で示される演算により、修正係数行列
0006905296.tif 000087
に変換し、前記（３）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
【数６】
0006905296.tif 000088
【数７】
0006905296.tif 000089
として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
【数８】
0006905296.tif 000090
で示される出力方程式を出力し、
前記状態方程式計算手段は、前記（６）式と（７）式とで示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｚ（ｔ）を計算し、更に、前記（８）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項７】
前記方程式設定手段は、
0006905296.tif 000091
を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
【数９】
0006905296.tif 000092
を設定し、
前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
【数１０】
0006905296.tif 000093
で示される演算により、固有値λ１、・・・、λｓ、・・・、λｎ（１≦ｓ≦ｎ）を前記対角要素とする前記固有値包含対角行列である対角行列Λに変換し、
前記固有値変更手段は、前記（１０）式で示される前記対角行列Λ内の各前記固有値λｓ＝σｓ＋ωｓｉ（１≦ｓ≦ｎ）が所定の条件を満たす場合に、前記固有値λｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
0006905296.tif 000094
を出力し、
前記方程式変換手段は、前記係数行列Ａを、
【数１１】
0006905296.tif 000095
で示される演算により、修正係数行列
0006905296.tif 000096
に変換し、前記係数行列Ｂを、
【数１２】
0006905296.tif 000097
で示される演算により、修正係数行列
0006905296.tif 000098
に変換し、前記（９）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
【数１３】
0006905296.tif 000099
として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
【数１４】
0006905296.tif 000100
で示される出力方程式を出力し、
前記状態方程式計算手段は、前記（１３）式で示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｘ（ｔ）を計算し、更に、前記（１４）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項８】
前記方程式設定手段は、
0006905296.tif 000101
を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
【数１５】
0006905296.tif 000102
を設定し、
前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
【数１６】
0006905296.tif 000103
【数１７】
0006905296.tif 000104
で示される演算により、固有値λ１、・・・、λｓ、・・・、λｍ（１≦ｓ≦ｍ）を前記対角要素に含む前記固有値包含対角行列であるジョルダン標準形行列Ｊに変換し、
前記固有値変更手段は、前記（１６）式及び前記（１７）式で示される前記ジョルダン標準形行列Ｊ内の各前記固有値λｓ＝σｓ＋ωｓｉ（１≦ｓ≦ｍ）が所定の条件を満たす場合に、前記固有値λｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
0006905296.tif 000105
を出力し、
前記方程式変換手段は、前記係数行列Ｂを、
【数１８】
0006905296.tif 000106
で示される演算により、修正係数行列
0006905296.tif 000107
に変換し、前記（１５）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
【数１９】
0006905296.tif 000108
【数２０】
0006905296.tif 000109
として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
【数２１】
0006905296.tif 000110
で示される出力方程式を出力し、
前記状態方程式計算手段は、前記（１９）式と（２０）式とで示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｚ（ｔ）を計算し、更に、前記（２１）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項９】
前記方程式設定手段は、
0006905296.tif 000111
を前記変数ベクトル項の１階導関数、Ｘ（ｔ）を前記変数ベクトル項、Ｕ（ｔ）を前記入力ベクトル項、Ａを前記変数ベクトル項の係数行列、Ｂを前記入力ベクトル項の係数行列として、前記線形状態方程式として、
【数２２】
0006905296.tif 000112
を設定し、
前記固有値計算手段は、前記係数行列Ａのサイズをｎを自然数としてｎ行ｎ列とし、正則な行列Ｐを用いて、前記係数行列Ａを、
【数２３】
0006905296.tif 000113
【数２４】
0006905296.tif 000114
で示される演算により、固有値λ１、・・・、λｓ、・・・、λｍ（１≦ｓ≦ｍ）を前記対角要素に含む前記固有値包含対角行列であるジョルダン標準形行列Ｊに変換し、
前記固有値変更手段は、前記（２３）式及び前記（２４）式で示される前記ジョルダン標準形行列Ｊ内の各前記固有値λｓ＝σｓ＋ωｓｉ（１≦ｓ≦ｍ）が所定の条件を満たす場合に、前記固有値λｓを所定の定数に変更することにより、前記修正固有値包含対角行列である修正対角行列として
0006905296.tif 000115
を出力し、
前記方程式変換手段は、前記係数行列Ａを、
【数２５】
0006905296.tif 000116
で示される演算により、修正係数行列
0006905296.tif 000117
に変換し、前記係数行列Ｂを、
【数２６】
0006905296.tif 000118
で示される演算により、修正係数行列
0006905296.tif 000119
に変換し、前記（２２）式で示される前記線形状態方程式に対応する前記修正線形状態方程式を、
【数２７】
0006905296.tif 000120
として出力し、更に、Ｃ及びＤを係数行列として、前記制御装置への前記出力信号のベクトルＹ（ｔ）を出力するための、
【数２８】
0006905296.tif 000121
で示される出力方程式を出力し、
前記状態方程式計算手段は、前記（２７）式で示される前記出力状態方程式を前記差分方程式の形式に変換し、前記差分方程式に対する前記解法演算を実行することにより解Ｘ（ｔ）を計算し、更に、前記（２８）式で示される前記出力方程式を他の差分方程式に変換し、前記他の差分方程式に対する他の解法演算を実行することにより前記出力信号のベクトルＹ（ｔ）を計算して前記制御装置に出力する、
請求項１乃至５の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項１０】
前記状態方程式計算手段は、Ｎを任意の自然数として、
【数２９】
0006905296.tif 000122
【数３０】
0006905296.tif 000123
で示される前記差分方程式に対する、Ｔａｙｌｏｒ展開法又はこれと等価な計算式となるＲｕｎｇｅ−Ｋｕｔｔａ法による解法演算を実行し、
或る計算ステップｈに対して、前記固有値λｓが、
【数３１】
0006905296.tif 000124
で示される条件を充足することを判別するゼロ点検出手段を更に備え、
前記固有値変更手段は、前記ゼロ点検出手段が前記（３１）式で示される条件が充足されたことを判別した場合に、前記固有値λｓを他の固有値に変更する、
請求項１、２、３、４、５、７、又は９の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項１１】
前記状態方程式計算手段は、Ｎを任意の自然数として、
【数３２】
0006905296.tif 000125
【数３３】
0006905296.tif 000126
で示される前記差分方程式に対する、Ｔａｙｌｏｒ展開法又はこれと等価な計算式となるＲｕｎｇｅ−Ｋｕｔｔａ法による解法演算を実行し、
或る計算ステップｈに対して、前記固有値λｓが、
【数３４】
0006905296.tif 000127
で示される条件を充足することを判別するゼロ点検出手段を更に備え、
前記方程式変換手段は、前記ゼロ点検出手段が前記（３５）式で示される条件が充足されたことを判別した場合に、前記Ｎの値を変更する、
請求項１、２、３、４、５、７、又は９の何れか一項に記載の制御対象の動作を模擬する計算装置。
【請求項１２】
制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置に、
制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定処理と、
前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算処理と、
前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更処理と、
前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換処理と、
前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算処理と、
を実行させる、計算装置において状態方程式を解く演算の実行方法。
【請求項１３】
制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、前記制御対象の動作を模擬する計算装置に、
制御装置への出力を表す変数ベクトル項と前記制御装置からの入力を表す入力ベクトル項を含み、前記変数ベクトル項及び前記入力ベクトル項に夫々乗算される各係数行列の要素が一定時間定数とみなせる線形時間不変システムに前記制御対象をモデル化した線形状態方程式を、前記状態方程式として設定する方程式設定処理と、
前記変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する固有値計算処理と、
前記固有値包含対角行列内の各前記固有値が所定の条件を満たす場合に、前記固有値の実部又は虚数部を所定の定数に変更することにより、前記固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更処理と、
前記線形状態方程式を前記修正固有値包含対角行列に基づく形式に変換することにより、前記線形状態方程式に対応する修正線形状態方程式を出力する方程式変換処理と、
前記修正線形状態方程式を差分方程式の形式に変換し、前記差分方程式に対する解法演算を実行することにより、前記制御装置への出力信号を生成する状態方程式計算処理と、
を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、制御対象の動作を数学的にモデル化した状態方程式を解く演算を実行することにより制御対象の動作を模擬し、それにより制御装置の設計及び検証を支援する計算装置、そのような計算装置において状態方程式を解く演算の実行方法、及びプログラムに関する。
【背景技術】
【０００２】
制御対象を接続することにより、その制御対象に対してフィードバック制御を行う制御装置が知られている。その制御対象は、機械装置である場合もあり、或いは電子回路である場合もある。例えば、自動車には、ＥＣＵ（Ｅｌｅｃｔｒｏｎｉｃｓ  Ｃｏｎｔｒｏｌ  Ｕｎｉｔ：電子制御ユニット）と呼ばれる制御装置が、数十個搭載されている。夫々のＥＣＵは、例えばエンジン、トランスミッション、又はサスペンションなどの機械装置である制御対象に対して、フィードバック制御を行う。また、ＥＣＵは、例えばダッシュボードパネルに各種情報を表示する液晶計器類をはじめとする電子回路である制御対象に対して、フィードバック制御を行う。これらの制御装置は、制御対象を接続して初めて機能するという特性を持つ。従って、制御装置の開発においては、かつては、実際に制御対象を接続して実験／検証を行う必要があった。しかし、開発段階の制御装置を実際の制御対象に接続して検証を行った場合、制御装置による制御が不完全であることにより、制御対象が破壊されてしまうという事故が起こり得た。また、例えば月着陸船の制御装置のように、実際の使用環境での実験／検証が不可能な制御装置の開発においては、制御対象を実際に接続してその制御装置を開発することができなかった。或いは、開発段階において制御対象の実機が存在しない場合においても、制御対象を実際に接続してその制御装置を開発することができなかった。
【０００３】
このような問題を解決するために、従来、ＨＩＬＳ（Ｈａｒｄｗａｒｅ  Ｉｎ  ｔｈｅ  Ｌｏｏｐ  Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置が知られている（例えば特許文献１）。ＨＩＬＳでは、機械装置又は電子回路である制御対象が計算装置におけるソフトウェア処理によって模擬される。この場合、制御対象の動作は、例えば常微分方程式である状態方程式によって構成される数学モデルで表現される。また、計算装置は、電圧又は電流等の実際の信号を入出力する物理的入出力手段を装備する。そして、計算装置のこれらの物理的入出力手段に、制御装置の実際の物理的入出力手段が接続される。計算装置は、状態方程式を解く演算を実時間で実行することにより制御対象の動作を模擬する。そして、計算装置は実際の信号を制御装置に対して実時間で入出力することにより、ハードウェアである制御装置を実際に動作させてテストする。即ち、ＨＩＬＳは、リアルタイムシミュレータとして動作する。ＨＩＬＳにより、制御対象の実機を用意することなく、ＥＣＵなどの制御装置の設計及び検証を支援することが可能となる。
【０００４】
更に、ＳＩＬＳ（Ｓｏｆｔｗａｒｅ  Ｉｎ  ｔｈｅ  Ｌｏｏｐ  Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置も知られている（例えば特許文献１）。ＳＩＬＳでは、制御装置自体、及びその物理的入出力手段も計算装置によって模擬される。そして、ＳＩＬＳは、制御装置と制御対象の全体の動作を、計算装置が実行するソフトウェア処理によって模擬する。
【０００５】
ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、一般的に常微分方程式である状態方程式として表現することができる。特に、変数ベクトルに第１の係数行列を乗算し、入力ベクトルに第２の係数行列を乗算した各乗算結果の和が、変数ベクトルの１階の導関数に等しい線形状態方程式は、制御対象の挙動を表す数学モデルとして有用であり、非線形の場合であっても上記形式に近似して用いられることが多い。ここで、「変数ベクトル」とは、常微分方程式の解をいう。
【０００６】
ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、線形状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実時間で入出力する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。
【０００７】
連続系の状態方程式を解くための計算方法としては、その状態方程式を離散信号形式で表現することにより差分方程式として計算する数値計算手法が良く知られている。この解法の代表例としては、Ｒｕｎｇｅ−Ｋｕｔｔａ法、及びＴａｙｌｏｒ展開法などが知られている。
【０００８】
ところで、線形状態方程式における係数行列の要素が或る一定時間において不変な定数であれば、その方程式は、所謂ＬＴＩ（Ｌｉｎｅａｒ  Ｔｉｍｅ  Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルは、ＬＴＩシステムを仮定しても十分であることが多い。このようなＬＴＩシステムの状態方程式においては、係数行列を対角要素が互いに異なる固有値である対角行列に変換することができる。
【０００９】
このような対角行列に含まれる固有値は、一般に複素数である。そして、複数の固有値の大きさが著しく異なる状態方程式は、解が周波数の比較的高い成分と低い成分を含む、いわゆるＳｔｉｆｆ状態方程式と呼ばれる（例えば特許文献２）。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特許第５６９２７３９号公報
【特許文献２】特開２００５−０２５６５１号公報
【発明の概要】
【発明が解決しようとする課題】
【００１１】
上述のようなＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、前述したＲｕｎｇｅ−Ｋｕｔｔａ法又はＴａｙｌｏｒ展開法などによる差分方程式に変換して演算する場合について考える。このケースでは、差分方程式における離散時間の計算ステップと対角行列の固有値が適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（−∞）に発散することが知られている。特に、前述したＳｔｉｆｆ状態方程式を解く場合に、絶対値の大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくしなければならない。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に対しては、その計算ステップは過剰に小さくなりすぎる。この結果、計算ステップを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。
【００１２】
このような問題を解決するために、更に従来、複素平面のいわゆる左半平面全体を計算安定領域とする、後退Ｅｕｌｅｒ法又は台形法（Ｔｕｓｔｉｎ法）と呼ばれる解法が知られている。しかしながら、これらの解法ではＬＴＩシステムの場合であっても逆行列を計算する必要がある。一般に逆行列の計算負荷は、行列が“ｎ×ｎ”要素のサイズを有する場合、ｎの３乗に比例すると言われている。従って、ｎが大きなシステムの場合には、計算負荷が極めて高くなる。このため、Ｅｕｌｅｒ法や台形法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。
【００１３】
更に、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次、台形法は２次の近似精度となる。このため、固有値λの絶対値が大きい場合は、これらの解法は、安定動作範囲が広いため発散しないで計算できるものの、計算精度が著しく落ちるというあまり知られていない事実がある。
【００１４】
以上の考察からわかるように、計算装置において絶対値が大きな固有値を有する方程式を適切に解く演算を実行することは、従来困難であった。
【００１５】
そこで、本発明は、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とすることを目的とする。
【課題を解決するための手段】
【００１６】
態様の一例は、制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行することにより、制御対象の動作を模擬する計算装置であって、制御装置への出力を表す変数ベクトル項Ｙと制御装置からの入力を表す入力ベクトル項Ｕを含み、変数ベクトル項Ｘ及び入力ベクトル項Ｕに夫々乗算される各係数行列Ａ，Ｂの要素が一定時間において不変な定数とみなせる線形時間不変システムに制御対象をモデル化した線形状態方程式を、状態方程式として設定する方程式設定手段と、変数ベクトル項Ｘの係数行列Ａを、対角要素が固有値である対角行列Λ又はジョルダン標準形行列Ｊの何れかの行列である固有値包含対角行列に変換する固有値計算手段と、固有値包含対角行列内の各固有値が所定の条件を満たす場合に、固有値の実部又は虚数部の一方又は両方を所定の定数に変更することにより、固有値包含対角行列に対応する修正固有値包含対角行列を出力する固有値変更手段と、線形状態方程式を修正固有値包含対角行列に基づく形式に変換することにより、線形状態方程式に対応する修正線形状態方程式を出力する方程式変換手段と、修正線形状態方程式を差分方程式の形式に変換し、差分方程式に対する解法演算を実行することにより、制御装置への出力信号を生成する状態方程式計算手段と、を備える。
【発明の効果】
【００１７】
本発明によれば、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、ＨＩＬＳまたはＳＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。
【図面の簡単な説明】
【００１８】
【図１】Ｒｕｎｇｅ−Ｋｕｔｔａ法において安定に計算できる固有値λと計算ステップhの積ｈλの複素平面領域を説明する図である。
【図２】Ｔａｙｌｏｒ展開法において安定に計算できる固有値λと計算ステップｈの積ｈλの領域を説明する図である。
【図３】本発明の第１の実施形態であるＨＩＬＳ計算装置の実施形態を示すブロック図である。
【図４】制御対象の機構例を示す図である。
【図５】制御対象のブロック線図の例を示す図である。
【図６】固有値実部変更手段及び固有値虚数部変更手段の動作例を示すフローチャートである。
【図７】本発明の第２の実施形態であるＳＩＬＳ計算装置の実施形態を示すブロック図である。
【図８】固有値λ＝±πｉを異なる計算手法で計算した時系列結果を説明する図である。
【図９】固有値λ＝±π／２ｉを異なる計算手法で計算した時系列結果を説明する図である。
【図１０】指数関数の時間応答を異なる計算手法で計算した時系列結果を説明する図である。
【図１１】複素数の固有値を異なる計算手法で計算した時系列結果を説明する図である。
【図１２】本発明の第３の実施形態であるＨＩＬＳ計算装置の実施形態を示すブロック図である。
【図１３】第３の実施形態における固有値λの変更方法を説明する図である。
【図１４】固有値変更手段の動作例を示すフローチャートである。
【図１５】４次のＴａｙｌｏｒ展開法におけるゼロ点の複素平面上での位置を説明する図である。
【図１６】固有値が−０．２７±２．５０ｉの場合において４次のＴａｙｌｏｒ展開法で計算した時の時系列計算結果と、第４の実施形態による改善例を説明する図である。
【発明を実施するための形態】
【００１９】
以下、本発明を実施するための形態（以下「本実施形態」と記載）について図面を参照しながら詳細に説明する。以下の説明においては、本実施形態について詳細に説明する前に、まず状態方程式に関する一般的な技術背景について説明する。そして、その一般技術の問題点を具体的に明らかにした上で、本実施形態について詳細に説明する。
【００２０】
前述したように、ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、任意数の要素を持つ状態変数ベクトルＸ、時間を表す変数ｔ、及びベクトル関数Ｆ（ｔ，Ｘ）を用いて、下記（１）式の常微分方程式として表現される。
【数１】
0006905296.tif 000002
【００２１】
上記（１）式の常微分方程式は、公知の技術である動作点による区分線形化手法などにより、Ａ、Ｂ行列と、外力を表す入力ベクトルＵ（ｔ）とを用いた、線形状態方程式と呼ばれる下記（２）式の形式で表現することができる。
【数２】
0006905296.tif 000003
【００２２】
このように、線形状態方程式は、変数ベクトルＸ（ｔ）に係数行列Ａを乗算し、入力ベクトルＵ（ｔ）に係数ベクトルＢを乗算した各乗算結果ＡＸ（ｔ）とＢＵ（ｔ）の和が、変数ベクトルの１階の導関数
0006905296.tif 000004
に等しい、という形を有する。
【００２３】
ＨＩＬＳまたはＳＩＬＳにおいて（２）式の形式の状態方程式を解く演算を効率的に実行することは重要である。特に、ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、（２）式の状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実現時間で供給する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。
【００２４】
上述の（２）式に示されるようなＡ、Ｂ行列で表現される連続系方程式を解くための計算方法としては、（２）式の形式の状態方程式を離散信号形式で表現することにより差分方程式として計算する数値計算手法が知られている。前述したように、この解法の代表例として従来、Ｒｕｎｇｅ−Ｋｕｔｔａ法、又はＴａｙｌｏｒ展開法などが知られている。
【００２５】
ところで、上記（２）式の各係数行列Ａ、Ｂは必ずしも定数行列である必要はない。しかしながら、各係数行列Ａ、Ｂの要素が一定時間において不変な定数であれば、（２）式は、所謂ＬＴＩ（Ｌｉｎｅａｒ  Ｔｉｍｅ  Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの線形状態方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルはＬＴＩシステムを仮定しても十分であることが多い。係数行列Ａ、Ｂが定数であるＬＴＩシステムでは、（２）式は下記（３）、（４）式に示されるように固有値分解することができる。
【数３】
0006905296.tif 000005
【数４】
0006905296.tif 000006
【００２６】
上記（４）式に示されるように、Λは対角要素が固有値である対角行列である。従って、（２）式は、スカラ連立方程式として扱うことができる。
【００２７】
上記（４）式の固有値λｓ（１≦ｓ≦ｎ）は、一般に複素数である。そして、前述したように、複数の固有値λｓの大きさが著しく異なる状態方程式は、いわゆるＳｔｉｆｆ状態方程式と呼ばれる。
【００２８】
なお、係数行列Ａを上記（４）式のように対角行列に分解ができない場合もある。その場合には、後述するように、係数行列Ａはジョルダン（Ｊｏｒｄａｎ）標準形行列に分解することができる。
【００２９】
上述の（３）、（４）式に示されるＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、前述したＲｕｎｇｅ−Ｋｕｔｔａ法又はＴａｙｌｏｒ展開法などの差分方程式に変換して演算する場合について考える。このケースでは、差分方程式における離散時間の計算ステップｈと、Ａ行列の固有値λの値とが適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（−∞）に発散することが知られている。逆に発散せず安定に状態方程式を解くことができる条件は、計算アルゴリズムごとに異なる。
【００３０】
図１は、Ｒｕｎｇｅ−Ｋｕｔｔａ法において安定に計算できる固有値λと計算ステップｈの積ｈλの複素平面領域を説明する図である。図１において、横軸は積ｈλの値の実部を示す実部軸、縦軸は積ｈλの値の虚数部を示す虚数部軸である。図１の安定領域曲線１０１、１０２、及び１０３は夫々、前述の（２）式の状態方程式を代表的な数値解法である１次、２次、及び４次のＲｕｎｇｅ−Ｋｕｔｔａ法で解いたときの、発散せずに計算できる安定境界において積ｈλの値をプロットしたものである。積ｈλの計算において、ｈは上記Ｒｕｎｇｅ−Ｋｕｔｔａ法が実行されるときの計算ステップであり、λは（２）式の係数行列Ａと正則な行列Ｐを用いて前述の（３）式によって計算される。
【００３１】
例えば、４次のＲｕｎｇｅ−Ｋｕｔｔａ法は、下記（５）式で示される。
【数５】
0006905296.tif 000007
【００３２】
図１で、各次数でのＲｕｎｇｅ−Ｋｕｔｔａ法の計算において、各次数に対応する安定領域曲線１０１、１０２、又は１０３の各内側に積ｈλが存在すれば、対応する次数のＲｕｎｇｅ−Ｋｕｔｔａ法により発散しない計算が可能である。
【００３３】
図１において、破線領域１０４は、Ｓｔｉｆｆな固有値に基づく積ｈλの存在領域である。また、破線領域１０５は、応答を支配する固有値に基づく積ｈλの存在領域である。前述したＳｔｉｆｆ状態方程式を解く場合に、絶対値が大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくして、積ｈλの存在領域１０４が安定領域曲線１０１、１０２、又は１０３の内側に入るようにする必要がある。計算ステップｈを小さくすれば、図１において、積ｈλの存在領域が領域１０４から破線矢印１０６で示される方向に移動して安定領域曲線１０１、１０２、又は１０３の内側に入るため、発散しない計算が可能である。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に基づく積ｈλの存在領域１０５付近では、その計算ステップｈは過剰に小さくなりすぎる。この結果、計算ステップｈを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。
【００３４】
前述した（２）式が前述したＬＴＩシステムの状態方程式となっており、かつ入力が一定値ｕであるとみなせる場合、行列指数関数
0006905296.tif 000008
を用いて、下記（６）式の一般解を得ることができる。
【数６】
0006905296.tif 000009
【００３５】
上記（６）式は、下記（７）式及び（８）式において、Ｎ＝∞（無限大）の場合に相当する。
【数７】
0006905296.tif 000010
【数８】
0006905296.tif 000011
【００３６】
しかし、Ｎ＝∞（無限大）の場合には、計算が困難である。このため、（７）式及び（８）式において有限のＮを用いると共に、計算誤差を少なくするために微小時間＝ｈの計算ステップで（７）式及び（８）式の差分方程式を計算する手法がある。この場合、ＸｎはＸ（ｈ×ｎ）の近似値であり、ｕは［ｔ×ｈ，（ｔ＋１）×ｈ］の間一定であるとし、ｕｎとする。上記（７）式及び（８）式の演算法を、本明細書ではＴａｙｌｏｒ展開法と呼ぶ。
【００３７】
図２は、Ｔａｙｌｏｒ展開法において安定に計算できる固有値λと計算ステップｈの積ｈλの複素平面領域を説明する図である。図２において、図１の場合と同様に、横軸は積ｈλの値の実部を示す実部軸、縦軸は積ｈλの値の虚数部を示す虚数部軸である。図２の安定領域曲線２０１、２０２、２０３、２０４、２０５、及び２０６は夫々、前述の（７）式及び（８）式においてN=１、２、４、５、８及び１０とした場合に計算したときの、発散せずに計算できる安定境界において積ｈλの値をプロットしたものである。積ｈλの計算において、ｈは上記（７）式及び（８）式の状態方程式の解法演算が実行されるときの計算ステップｈに等しく、λは（７）式及び（８）式の係数行列Ａと正則な行列Ｐを用いて前述の（３）式によって計算される。
【００３８】
なお、図２において、積ｈλの虚数部の値が−πｉから＋πｉ（ｉは虚数単位）までの間のグレー色に塗りつぶされた範囲２０７は、ナイキスト周波数を超えない範囲であり、その外側の範囲はナイキスト周波数を超える範囲である。積ｈλとナイキスト周波数との関係については、本実施形態の説明において後述する。
【００３９】
また、図２において、図１の場合と同様に、各次数においてでのＴａｙｌｏｒ展開法の計算において、各次数に対応する安定領域曲線２０１、２０２、２０３、２０４、２０５、又は２０６の各内側に積ｈλが存在すれば、発散しない計算が可能である。なお、Ｔａｙｌｏｒ展開法の次数とは前述した（８）式におけるＮの値を指すものである。
【００４０】
なお、（７）式及び（８）式が前述したＬＴＩシステムの状態方程式となっている場合には、前述した（５）式のＲｕｎｇｅ−Ｋｕｔｔａ法は、式を整理すれば４次のＴａｙｌｏｒ展開法となる。４次までのＲｕｎｇｅ−Ｋｕｔｔａ法の安定領域（図１の１０１、１０２、及び１０３）は、Ｔａｙｌｏｒ展開法の安定領域（図２の２０１、２０２、及び２０３）と同一であることが、一般に知られている。
【００４１】
以上説明したＲｕｎｇｅ−Ｋｕｔｔａ法及びＴａｙｌｏｒ展開法は、広く使われる数値計算手法である。しかし、これらの数値計算法は、安定領域に限りがある。そして、Ｓｔｉｆｆ状態方程式の場合、前述したように、不安定になり計算が発散するか、過剰に小さな計算ステップｈを用いなければならないという問題がある。
【００４２】
このような問題を解決するために、複素平面のいわゆる左半平面全体を計算安定領域とする、下記（９）式で示される後退Ｅｕｌｅｒ法又は下記（１０）式で示される台形法（Ｔｕｓｔｉｎ法）と呼ばれる解法が、一般的に用いられている。
【数９】
0006905296.tif 000012
【数１０】
0006905296.tif 000013
【００４３】
しかしながら、（９）式又は（１０）式に示される解法では逆行列を計算する必要がある。一般に逆行列の計算負荷は、行列が“ｎ×ｎ”要素のサイズを有する場合、ｎの３乗に比例すると言われている。従って、ｎが大きなシステムの場合には、計算負荷が極めて高くなる。このため、Ｅｕｌｅｒ法や台形法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。
【００４４】
更に、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次、台形法は２次の近似精度となる。このため、固有値λの絶対値が大きい場合は、これらの解法は、安定動作範囲が広いため発散しないで計算できるものの、計算精度が著しく落ちるというあまり知られていない事実がある。
【００４５】
以上の考察からわかるように、計算装置において絶対値が大きな固有値λを有する方程式を適切に解く演算を実行することは、一般的に困難である。
【００４６】
そこで、本実施形態は、制御対象の動作を模擬する計算装置において、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とすることを目的とする。本実施形態の詳細について、以下に説明する。
【００４７】
本実施形態の動作原理について説明する。まず、一般的な状態方程式である前述の（２）式において、変数ベクトル項Ｘ（ｔ）の係数行列Ａを、対角要素が固有値である対角行列、又はジョルダン標準形行列に変換する。本実施形態において、これらの行列を総称して、固有値包含対角行列と呼ぶことにする。
【００４８】
例えば、（２）式の係数行列Ａ、Ｂが定数であるＬＴＩシステムにおいて、下記（１１）式（前述した（３）式の再掲）に示される変換演算が実行される。
【００４９】
【数１１】
0006905296.tif 000014
【００５０】
この結果、（２）式の固有値を対角要素に持つ下記（１１）式（前述の（４）式の再掲）に示される対角行列Λが、上述の固有値包含対角行列として算出される。
【数１２】
0006905296.tif 000015
【００５１】
或いは、（２）式の係数行列Ａのサイズをｎ行ｎ列（ｎは自然数）とし、正則な行列Ｐを用いて、係数行列Ａに対して、下記（１３）式及び（１４）式で示される演算が実行される。この結果、固有値λ１、・・・、λｓ、・・・、λｍ（１≦ｓ≦ｍ）を対角要素に持つ（１３）式で示されるジョルダン細胞行列と呼ばれる行列を、更に対角要素として持つ（１４）式で示されるジョルダン標準形行列Ｊが算出される。結果的に、ジョルダン標準形行列Ｊは、その対角要素ｓ（１≦ｓ≦ｍ）に固有値λｓを含む、固有値包含対角行列である。
【数１３】
0006905296.tif 000016
【数１４】
0006905296.tif 000017
【００５２】
そして、本実施形態では、上記固有値包含対角行列内の各固有値が所定の条件を満たす場合に、その固有値の実部又は虚数部の一方又は両方を所定の定数に変更する。
【００５３】
例えば、対角行列Λ内の固有値λｓ（１≦ｓ≦ｎ）、又はジョルダン標準形行列Ｊ内の固有値λｓ（１≦ｓ≦ｍ）の実部σｓの絶対値｜σｓ｜が所定の実部閾値ａより大きい場合に、上記固有値λｓの実部σｓを所定の実数定数ｂに変更する。より具体的には、任意の数値計算手法、例えば前述の（７）式及び（８）式のＴａｙｌｏｒ展開法の演算に着目する。例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば実数“−３”の絶対値“３”を、上記実部閾値ａとする。そして、固有値λｓの実部σｓの絶対値｜σｓ｜が上記実部閾値ａ＝３より大きい場合に、固有値λｓの実部σｓを所定の実数定数ｂに変更するものである。
【００５４】
また例えば、対角行列Λ内の固有値λｓ（１≦ｓ≦ｎ）、又はジョルダン標準形行列Ｊ内の固有値λｓ（１≦ｓ≦ｍ）の虚数部ωｓｉの絶対値｜ωｓ｜が所定の虚数部閾値ｃより大きい場合に、上記固有値λｓの虚数部ωｓｉを所定の虚数定数ｄｉに変更する。例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば虚数“πｉ”（ｉ：虚数単位）の絶対値“π”を、上記虚数部閾値ｃとする。そして、固有値λｓの虚数部ωｓｉの絶対値｜ωｓ｜が上記虚数部閾値ｃ＝πより大きい場合に、固有値λｓの虚数部ωｓｉを所定の虚数定数ｄｉ（ｄは正の実数）に変更するものである。また例えば５次のＴａｙｌｏｒ展開法の安定領域内（図２の安定領域曲線２０５の範囲内）に存在する例えば虚数“−πｉ”の絶対値“π”を、上記虚数部閾値ｃ′とする。このとき、固有値λｓの虚数部ωｓｉの絶対値｜ωｓ｜が上記虚数部閾値ｃ′＝πより大きい場合に、固有値λｓの虚数部ωｓｉを所定の虚数定数−ｄｉ（ｄは正の実数）に変更するものである。
【００５５】
従って、例えば（１２）式の対角行列Λ内の各対角要素λ２及びλ３において、その実部の絶対値｜σ２｜及び｜σ３｜が実部閾値ａ＝３より大きく、かつその虚数部の絶対値｜ω２｜及び｜ω３｜が虚数部閾値ｃ＝π又はｃ′＝πより大きければ、対角行列Λが下記（１５）式に例示されるように変更されるものである。この行列を、修正対角行列と呼ぶ。修正対角行列は、固有値包含対角行列を変更した修正固有値包含対角行列に対応する。
【数１５】
0006905296.tif 000018
【００５６】
その後、（２）式に例示される状態方程式が固有値包含対角行列内の固有値を変更して得られる修正固有値包含対角行列に基づく形式に変換される。そして、変換後の状態方程式に基づいて、例えば１０次のＴａｙｌｏｒ展開法により得られる差分方程式に対する解法演算が実行される。
【００５７】
以上の動作原理によれば、（２）式に例示される状態方程式を常に積ｈλに関する安定領域内で解くことが可能となる。このため、状態方程式がたとえＳｔｉｆｆ状態方程式であったとしても、上記安定性が失われることはない。そして、状態方程式が変更された後は、在来のＲｕｎｇｅ−Ｋｕｔｔａ法やＴａｙｌｏｒ展開法により得られる差分方程式に対する解法演算が実行される。この結果、本出願によれば、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、ＨＩＬＳおよびＳＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。
【００５８】
上記動作原理に基づく本発明のいくつかの実施形態の詳細について、以下に説明する。図３は、本発明の第１の実施形態である、ＨＩＬＳ計算装置３００の実施形態を示すブロック図である。ＨＩＬＳ計算装置３００は、物理的な信号入出力手段を有する、いわゆるＨＩＬＳの実施形態である。
【００５９】
ＨＩＬＳ計算装置３００は、特には図示しないが、所謂組込みコンピュータのハードウェアとして、ＣＰＵ（中央演算処理装置）、ＲＯＭ（リードオンリーメモリ）、ＲＡＭ（ランダムアクセスメモリ）、フラッシュメモリ装着装置などを有する。図３に示されるＨＩＬＳ計算装置３００において、通信手段３０１は所謂イーサネットあるいはＣＡＮ通信、またはこれに準ずるものである。固有値計算手段３０２、固有値実部変更手段３０３、固有値虚数部変更手段３０４、行列算出手段３０５、ベクトル算出手段３０８、及び状態方程式計算手段３０９は、上記ＣＰＵが上記フラッシュメモリ装着装置に装着されたメモリカードに記憶されたＨＩＬＳ制御プログラムを上記ＲＡＭに読み出して実行する機能ブロック群である。時計装置３１０は、ＨＩＬＳ計算装置３００のコンピュータが内蔵するタイマ装置のハードウェアである。デジタル信号入力手段３０６、電圧入力手段３０７、電圧出力手段３１１は、上記入出力インタフェース装置及び上記ＲＡＭ又は上記外部記憶装置によって構成される。通信手段３０１には、ローカルエリアネットワークやインターネット等のネットワークまたは専用通信線を介してホストコンピュータ３２０が通信接続される。
【００６０】
デジタル信号入力手段３０６及び電圧入力手段３０７は、計算装置側物理的入力手段として動作し、特には図示しない制御装置の制御装置側物理的出力手段のデジタル信号出力手段及び電圧出力手段に接続される。
【００６１】
デジタル信号入力手段３０６は、制御装置からの複数の信号線の各電圧を監視し、夫々所定の閾値を超えた場合に“１”、超えない場合に“０”の各デジタル入力信号をベクトル算出手段３０８に出力する。
【００６２】
電圧入力手段３０７は、複数のアナログ／デジタル変換器を内蔵し、制御装置からの複数のアナログ入力電圧値をデジタル入力電圧値に変換し、各デジタル入力電圧値をベクトル算出手段３０８に出力する。
【００６３】
ベクトル算出手段３０８は、上記複数のデジタル入力信号の“１”又は“０”値及び複数のデジタル入力電圧値をベクトル要素とする入力ベクトルＵ（前述した（２）式のＵ（ｔ）に対応）を算出する。なお、入力ベクトルＵの第１要素としては値“１”を設定する。この第１要素は、状態方程式が定数項（前述した（２）式のＢＵ（ｔ）＝定数）を必要とした場合に、係数行列Ｂの値により任意の定数を生成するために使用する。
【００６４】
一方、ホストコンピュータ３２０は、一般的なパーソナルコンピュータであってよく、ＨＩＬＳの数学モデルを生成するためのシミュレーションモデリングツールのプログラムを実行することができる。ユーザは、例えば図４に例示されるような機械システムを制御対象とし、その数学モデルを生成したい場合には、上記シミュレーションモデリングツールを実行することにより、例えば図５に例示されるようなブロック線図を、ホストコンピュータ３２０のキーボード、マウス、及びディスプレイを使って作成する。
【００６５】
ホストコンピュータ３２０には更に、例えばシミュレーションモデリングツールの機能として、状態方程式生成ソフトウェアが搭載されている。この状態方程式生成ソフトウェアは、ユーザが図５に例示されるように作成したブロック線図に記載されている積分器の出力値を要素とする、例えば下記（１６）式として例示されるような状態変数ベクトルを生成する。
【数１６】
0006905296.tif 000019
【００６６】
状態方程式生成ソフトウェアは更に、上記ブロック線図を解析し、積分器に与えられる入力に基づいて、係数行列Ａ及びＢを生成する。例えば、図５に例示されるブロック線図におけるｖ１の出力である第１積分器５０１に接続されている値は、第２積分器５０２及び第３積分器５０３に夫々所定の係数を乗じてｍ１で除算したものである。従って、状態方程式生成ソフトウェアは、下記（１７）として例示される係数行列Ａ及びＢを生成する。
【数１７】
0006905296.tif 000020
【００６７】
上記のようにして算出される係数行列Ａ及びＢにより、下記（１８）式の状態方程式が構成されることになる。
【数１８】
0006905296.tif 000021
【００６８】
ここで、Ｆは、制御装置から与えられる入力値である。従って、（１８）式の状態方程式は、前述した（２）式の状態方程式に対応する。
【００６９】
更に、状態方程式生成ソフトウェアは、ユーザが指定することにより、制御装置に接続したい出力値をＹベクトルとして設定することができ、このＹベクトルを生成するための下記（１９）式として示される出力方程式を構成する係数行列Ｃ及びＤを生成する。
【数１９】
0006905296.tif 000022
【００７０】
なお、図1で示したＳｔｉｆｆな固有値の例（図１の領域１０４を参照）は、図４に例示した機構に対して、以下の諸元を与えたものである。
ｍ１＝３００  ｍ２＝３００  ｋｓ＝１０  ｋ＝１０  ｃ＝１０００
このように、Ｓｔｉｆｆ状態方程式は極めて一般的に存在する。特に、ＳＩＬＳやＨＩＬＳの用途である機械装置や電子回路の模擬においては、いままで避けて通れない難題であった。本実施形態は、これを解決する手段を提供するものである。
【００７１】
以上説明したようにしてホストコンピュータ３２０で生成された係数行列Ａ，Ｂ，Ｃ，及びＤのデータが、ホストコンピュータ３２０からＨＩＬＳ計算装置３００に送られ、図３の通信手段３０１で受信され、記憶される。
【００７２】
図３の説明に戻り、固有値計算手段３０２は、通信手段３０１が受信した上記係数行列Ａを受け取り、この係数行列Ａに基づいて下記の手順により固有値を算出する。
【００７３】
手順１として、下記（２０）式として示される固有方程式を計算する。
【数２０】
0006905296.tif 000023
係数行列Ａのサイズがｎ行ｎ列の場合、固有値λはｎ次多項式となる。ここで、ベアストウ法として知られた手法を用いることによって、全ての固有値λ１、・・・、λｓ、・・・、λｎを算出する。これらの固有値λｓ（１≦ｓ≦ｎ）に重複がなければ、これらを対角要素とする対角行列Λを、前述した（１２）式の形で生成する。
【００７４】
次に、固有値λｓ（１≦ｓ≦ｎ）に重複がなければ、手順２として、固有値λｓに対する固有ベクトルｘｓ（１≦ｓ≦ｎ）を、下記（２１）式により算出する。
【数２１】
0006905296.tif 000024
この計算は、ガウス消去法として知られる連立方程式の解を得るアルゴリズムにより実行される。
【００７５】
固有値λｓ（１≦ｓ≦ｎ）に重複がなければ、更に手順３とし、手順２で算出した固有ベクトルｘｓ（１≦ｓ≦ｎ）を用いて、下記（２２）式として示される正則行列Ｐを構成する。
【数２２】
0006905296.tif 000025
これと共に、Ｐの逆行列Ｐ−１を、ガウス・ジョルダン法と呼ばれる手法で計算する。
【００７６】
固有値λｓ（１≦ｓ≦ｎ）に重複がなければ、手順１で算出した対角行列Λを、固有値実部変更手段３０３及び固有値虚数部変更手段３０４に出力する。また、手順３で得られた、正則行列Ｐ及びその逆行列Ｐ−１と、上記対角行列の逆行列Λ−１とを、後述するように係数行列Ｂを修正するために、行列算出手段３０５に出力する。
【００７７】
なお、手順１で算出された固有値λｓ（１≦ｓ≦ｎ）に重複がある場合には、固有値計算手段３０２は、対角行列Λの代わりに、前述した（１３）式及び（１４）式で示したジョルダン標準形行列Ｊを算出する。ジョルダン標準形行列Ｊの算出方法は従来知られているが、ここではその説明は省略する。
【００７８】
また、手順１で算出された固有値λｓ（１≦ｓ≦ｎ）に重複がある場合には、重複した固有値に微小値を加減算することにより重複を解除して対角行列Λを作成して、ＳＩＬＳやＨＩＬＳに用いる場合は実用上の問題はない。
【００７９】
次に、固有値実部変更手段３０３及び固有値虚数部変更手段３０４の動作について説明する。図６は、固有値変更手段として一体として動作する固有値実部変更手段３０３及び固有値虚数部変更手段３０４の動作例を示すフローチャートである。
【００８０】
まず、図３の固有値計算手段３０２から、対角行列Λを取得する（図６のステップＳ６０１）。
【００８１】
次に、ｎ個の固有値λｓ（１≦ｓ≦ｎ）について処理を繰り返すための変数ｓに値１をセットする（図６のステップＳ６０２）。
【００８２】
その後、変数ｓの値が上記ステップＳ６０２で値１にセットされた後、図６のステップＳ６０９で変数ｓの値が＋１ずつインクリメントされながら図６のステップＳ６１０で変数ｓの値が値ｎを超えたと判断するまで、図６のステップＳ６０３からステップＳ６０８までの一連の処理を繰り返し実行する。以下、これら一連の処理について順次説明する。
【００８３】
まず、ステップＳ６０１で取得した対角行列Λから、変数ｓの値が示す位置の１つの対角要素の固有値“λｓ＝σｓ＋ωｓｉ”を取得する（図６のステップＳ６０３）。ここで、σｓは固有値λｓの実部、ωｓｉはｉを虚数単位として固有値λｓの虚数部である。
【００８４】
次に、ステップＳ６０３で取得した固有値λｓの実部σｓの絶対値｜σｓ｜が所定の実部閾値ａより大きいか否かを判定する（図６のステップＳ６０４）。
【００８５】
ステップＳ６０４の判定がＹＥＳならば、固有値λｓの実部σｓを所定の実数定数ｂに変更する（図６のステップＳ６０５）。
【００８６】
ステップＳ６０４の判定がＮＯならば、固有値λｓの実部σｓはそのままとして変更しない。
【００８７】
以上のステップＳ６０４とＳ６０５の処理が、図３の固有値実部変更手段３０３の機能を実現する。
【００８８】
次に、ステップＳ６０３で取得した固有値λｓの虚数部ωｓｉの絶対値｜ωｓ｜が所定の虚数部閾値ｃより大きいか否かを判定する（図６のステップＳ６０６）。
【００８９】
ステップＳ６０６の判定がＹＥＳならば、固有値λｓの虚数部ωｓｉを所定の虚数定数ｄｉ（ｉは虚数単位）に変更する（図６のステップＳ６０７）。
【００９０】
ステップＳ６０６の判定がＮＯならば、固有値λｓの虚数部ωｓｉはそのままとして変更しない。
【００９１】
以上のステップＳ６０６とＳ６０７の処理が、図３の固有値虚数部変更手段３０４の機能を実現する。
【００９２】
その後、適宜変更後の“σｓ＋ωｓｉ”を、修正対角行列（修正固有値包含対角行列）
0006905296.tif 000026
において、変数ｓの現在値が示す位置の対角要素の固有値λｓとして設定する（図６のステップＳ６０８）。
【００９３】
ステップＳ６０８の処理の後、変数ｓの値を＋１インクリメントする（図６のステップＳ６０９）。
【００９４】
そして、インクリメント後の変数ｓの値が対角行列Λの次数ｎを超えたか否かを判定する（図６のステップＳ６１０）。
【００９５】
ステップＳ６１０の判定がＮＯならば、ステップＳ６０３の処理に戻って、新たな変数ｓの値が示す対角行列Λ上の対角要素位置の固有値λｓの変更処理を実行する。
【００９６】
ステップＳ６１０の判定がＹＥＳになると、図３の固有値実部変更手段３０３及び固有値虚数部変更手段３０４の今回の計算ステップにおける処理を終了し、修正対角行列
0006905296.tif 000027
を、後述する図３の行列算出手段３０５及び状態方程式計算手段３０９に出力する。
【００９７】
以上のようにして、本実施形態では、修正対角行列
0006905296.tif 000028
に基づいて実行される状態方程式が常に積ｈλに関する安定領域内で解かれるように、対角行列Λの各対角要素の固有値を例えば前述した（１５）式に例示される修正対角行列の各対角要素のように変更することが可能となる。
【００９８】
図３の説明に戻り、行列算出手段３０５は、係数行列Ｂに関する以下の修正処理を実行する。まず、前述した（３）式及び（４）式で示されるＬＴＩシステムの状態方程式を前提として、下記（２３）式に示される方程式を設定する。（２３）式の方程式においては、t=∞値（無限大）を、対角行列Λで計算した場合と、修正対角行列
0006905296.tif 000029
で計算した場合に結果が定常値でありかつ同一であるという条件が設定される。
【数２３】
0006905296.tif 000030
上記（２３）式より、下記（２４）式が算出される。
【数２４】
0006905296.tif 000031
行列算出手段３０５は、上記数（２４）式の演算を実行する。即ち、図３において、行列算出手段３０５は、前述した図６のフローチャートに基づいて一体で動作する固有値実部変更手段３０３及び固有値虚数部変更手段３０４から、修正対角行列
0006905296.tif 000032
を入力する。また、行列算出手段３０５は、固有値計算手段３０２から、前述した正則行列Ｐ、その逆行列Ｐ−１、及び対角行列の逆行列Λ−１を入力する。更に、行列算出手段３０５は、通信手段３０１から、係数行列Ｂを入力する。行列算出手段３０５は、これらの入力したデータに基づいて、（２４）式により、係数行列Ｂから修正係数行列
0006905296.tif 000033
を算出し、状態方程式計算手段３０９に出力する。
【００９９】
最後に、図３の状態方程式計算手段３０９は、前述した（３）式及び（４）式で示されるＬＴＩシステムの状態方程式を前提として、以下に説明する状態方程式の計算処理を実行する。
まず、状態方程式計算手段３０９は、所定の値又はホストコンピュータ３２０からＨＩＬＳ計算装置３００内の通信手段３０１に転送された初期値ベクトル
を用いて、前述した（３）式に基づく下記（２５）式により、初期状態変数ベクトルＺ（０）を算出する。
【数２５】
0006905296.tif 000034
以後、状態方程式計算手段３０９は、前述の（３）式に基づく下記（２６）式により、状態変数ベクトルＺ（ｔ）を順次算出する。
【数２６】
0006905296.tif 000035
【０１００】
状態方程式計算手段３０９は、（２６）式の計算を実際には例えば、前述した（７）式及び（８）式を修正して得られる下記（２７）式及び（２８）式による、例えば５次のＴａｙｌｏｒ展開法の差分方程式として実行する。
【数２７】
0006905296.tif 000036
【数２８】
0006905296.tif 000037
【０１０１】
状態方程式計算手段３０９は、（２７）式及び（２８）式の差分方程式の解法演算を、時計装置３１０が定める所定の時間間隔ｈ秒毎に離散時間ｎをインクリメントしながら実行する。上記（２７）式及び（２８）式の演算において、入力値ｕｎは、ベクトル算出手段３０８から入力する入力ベクトルＵ（ｔ）を、時間間隔ｈでサンプリングして得られるものである。
【０１０２】
なお、Ｚｎは一般に複素数であるが、このため、状態方程式計算手段３０９は、（２７）式及び（２８）式により得られる状態変数ベクトルＺｎに対して、下記（２９）式として示される出力方程式を適用することにより、実数である出力ベクトルＹｎを得て電圧出力手段３１１に出力する。
【数２９】
0006905296.tif 000038
【０１０３】
電圧出力手段３１１は、状態方程式計算手段３０９における（２９）式の演算により得られた出力ベクトルＹｎに相当する電圧値を出力する。Ｙｎは時間間隔ｈごとに計算されるため、電圧出力手段３１１の出力は、時間間隔ｈで新しい電圧値に更新される。図３には図示していないが、電圧出力手段３１１内のデジタル／アナログ変換器によりデジタル電圧信号から変換されて出力されたアナログ電圧値は、制御装置の制御装置側物理的入出力手段に入力される。制御装置は、このアナログ電圧値を、制御対象に設置されたセンサ装置からの信号と理解して、制御対象の挙動を判定し、所定の制御出力を出力する。この制御出力は、ＨＩＬＳ計算装置３００のデジタル信号入力手段３０６又は電圧入力手段３０７に入力される。このようにして、制御装置と制御対象であるＨＩＬＳ計算装置３００との間に、フィードバックのクローズドループ（閉ループ）が形成される。
【０１０４】
図７は、図３とは異なる本発明の第２の実施形態の構成を示す図であり、ＳＩＬＳ計算装置７００の実施形態を示すブロック図である。図７は、計算装置がＳＩＬＳを構成する例である。図７では、図３の通信手段３０１の代わりに、行列設定手段７０２を有し、状態方程式において用いられる前述したＡ，Ｂ，Ｃ，Ｄ行列が設定される。図７において、図３の場合と同じ参照番号を付したブロックは、図３の場合と同じ機能を有する。
【０１０５】
また、ＳＩＬＳ計算装置７００は、制御装置制御則計算手段７０１を備えている。図３のＨＩＬＳ計算装置３００では、制御装置は、デジタル信号入力手段３０６、電圧入力手段３０７、又は電圧出力手段３１１を介してＨＩＬＳ計算装置３００の外部に接続されていた。これに対して、図７では、上記制御装置の動作を、ＳＩＬＳ計算装置７００内の制御装置制御則計算手段７０１によって模擬することができる。
【０１０６】
ただし、状態方程式計算手段３０９は計算ステップｈを用いて計算するが、制御装置制御則計算手段７０１の動作は必ずしも物理的にｈ秒毎に実行される必要はない（図７中では「ｈ′」と記載）。この点において、図７のＳＩＬＳ計算装置７００と図３のＨＩＬＳ計算装置３００は異なる。
【０１０７】
図７のＳＩＬＳ計算装置７００によれば、それ単体で、内蔵される制御装置制御則計算手段７０１によって制御装置の制御則の動作を検証することが可能である。電圧出力手段３１１の出力は、制御装置制御則計算手段７０１にフィードバックされると共に、ディスプレイ装置７０３に表示させて確認することができる。制御装置制御則計算手段７０１から状態方程式計算手段３０９への入力ベクトルＵ（ｔ）についても、ディスプレイ装置７０３で確認できるようにしてもよい。
【０１０８】
図３のＨＩＬＳ計算装置３００又は図７のＳＩＬＳ計算装置７００において、固有値実部変更手段３０３及び固有値虚数部変更手段３０４が固有値変更手段として動作する効果について、以下に説明する。
【０１０９】
前述した（２）式が固有値λを持つ場合、その解は前述した（６）式で与えられ、下記（３０）式となる。
【数３０】
0006905296.tif 000039
λが虚数部を持つ場合は、指数関数は下記（３１）式のオイラー公式より、時間的に振動する波形として現れる。
【数３１】
0006905296.tif 000040
【０１１０】
計算ステップｈ＝１の場合、これはサンプリング周波数＝１秒であるから、ナイキスト周波数は0.5Ｈｚ、周期は２秒であり、これ以上の周波数は表現できない。即ち、図２の領域２０７として示したように、ｈ＝１の計算ステップでは、周波数が0.5Ｈｚ＝πラジアン／秒以上である波形は、次数が高いＴａｙｌｏｒ展開法であれば安定領域に含まれる。この場合、発散しない計算が可能であるが、上記波形は、ｈ＝１のサンプリング周波数では表現されず、ナイキスト周波数で折り返して0.5Ｈｚより低い周波数、所謂エイリアス周波数となって現れる。
従って、固有値λの虚数部の絶対値がπラジアン／秒を超えた場合、この虚数部を例えば±πｉに変更して、その変更後の状態方程式を安定に計算できる５次のＴａｙｌｏｒ展開法で計算してもよい。
この場合、表現できない波形という意味を忠実に再現するならば、０ラジアン／秒としてもよい。
なお、Ｔａｙｌｏｒ展開法としては、前述した（７）式及び（８）式を用いることができる。そして、λ＝πｉの場合、例えば４次の項は、π４／２４≒４．０５９、５次の項はπ５／１２０≒２．５５となる。この結果は、Ｎ＝∞において得られる真値に対して、誤差が大きい。
これに対して、例えばＮ＝１０の場合、π１０／１０！＝０．０２５８となる。この場合、打ち切っても誤差が少ないと判断できる。このため、ｈλ＝πという大きな固有値をある程度の精度で計算するためには、この程度まで次数を上げる必要があり、５次では正しい値は得られない。
【０１１１】
図８は、状態方程式が、例えば下記（３２）式の方程式
【数３２】
0006905296.tif 000041
のように、固有値λ＝±πｉを持つ場合に、異なる計算手法で計算した時系列結果を説明する図である。図８は、ｈ＝１として、台形法、１０次のＴａｙｌｏｒ展開法、及び後退Ｅｕｌｅｒ法の計算結果を、夫々、８０１、８０２、及び８０３として、理論値８０４と共に示している。図８（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図８（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。図８から理解されるように、後退Ｅｕｌｅｒ法の計算結果８０３は、振動波形を表現できない。また、台形法の計算結果８０１は、時間軸対称という特性により２次の割には振動波形に強いが、やはり真値とはかけ離れている。
図８に関する考察より、例えばλ＝±πｉのような固有値は、確かにサンプル時間がｈ＝１の場合は丁度ナイキスト周波数であるため波形を表現し得る。しかし、このような固有値を持つ方程式を計算するためには、計算結果８０４として示される１０次のＴａｙｌｏｒ展開法のような、計算負荷が高い手法を採用せざるを得ない。
【０１１２】
そこで、本実施形態では、例えばサンプル時間h=1として、図３又は図７の固有値虚数部変更手段３０４において、前述した虚数部閾値ｃがπ／２（＝０．２５Ｈｚ）に設定される。そして、固有値の虚数部の絶対値が上記虚数部閾値ｃ＝π／２より大きなものについては、その固有値の虚数部が全て所定の虚数定数ｄｉ＝π／２ｉ又はｄｉ＝−π／２ｉに変更される。
図９は、固有値λ＝±π／２ｉを異なる計算手法で計算した時系列結果を説明する図である。図８の場合と同様に、図９（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図９（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。Ｔａｙｌｏｒ展開法の４次の項は、（π／２）４・２４＝０．２５４となる。このため、４次で打ち切られるＴａｙｌｏｒ展開法、又は４次のＲｕｎｇｅ−Ｋｕｔｔａ法を用いても、図９に示されるように、十分な計算精度を得ることができる。なお、図１及び図２に関して前述したように、前述した（２）式についてＬＴＩが仮定でき（３）式及び（４）式のように表現できるならば、４次以下の次数については、Ｔａｙｌｏｒ展開法とＲｕｎｇｅ−Ｋｕｔｔａ法は同一である。
この場合、図９のＲｕｎｇｅ−Ｋｕｔｔａ法の計算結果９０１は、本来の方程式が与える振動波形の周波数の理論値８０４とは明らかに異なる。しかし、従来用いられている代表的な手法である台形法の計算結果８０１であっても、本来の周波数の波形を計算することはできない。このため、少なくとも振動現象を表現できるという点では大差はなく、固有値λを±π／２ｉに変換して、Ｒｕｎｇｅ−Ｋｕｔｔａ法により計算した手法は計算効率に勝る。
【０１１３】
以上のように、そもそもにおいて例えばサンプル時間h=1の場合は固有値の虚数部の絶対値がπより大きければ、この固有値による応答は表現することができない。また、固有値の虚数部の絶対値がπ以下の値であっても、１０次のＴａｙｌｏｒ展開法のような特殊な手法を除き、実用になる数値計算法が存在せず、正確な波形再現はできない。そこで、本実施形態では、正確な計算が可能な位置まで固有値を変更することで、結果的により適切な計算結果を得ることがものである。
【０１１４】
次に固有値の実部の影響について説明する。図１０は、指数関数の時間応答を異なる計算手法で計算した結果を説明する図である。横軸は時間（秒）、縦軸はｘの計算結果である。図１０は、指数関数が負の実部を持つ場合であって、ｔ＝０での初期値を１とした場合の時間応答である。
ｅλｔの最初の１ステップ、即ちｔ＝ｈの値は、
λ＝−１  ｈ＝１：ｅ−ｔ  ＝０．３６８
λ＝−２  ｈ＝１：ｅ−２ｔ＝０．１３５
λ＝−３  ｈ＝１：ｅ−３ｔ＝０．０５０
λ＝−４  ｈ＝１：ｅ−４ｔ＝０．０１８３
である。従って、例えば固有値がλ＝−３より小さければ、その固有値の時間応答特性がλ＝−３のときの時間応答特性と同じである、として実用上は問題はない。２ステップ以降（ｔ＝２ｈ，３ｈ，４ｈ）では、時間応答特性は更に０に近似する。従って、固有値がλ＝−３より小さい場合と固有値がλ＝−３の場合とで、それらの差異は全くなくなる。
なお、“λ＝−３”のような原点から離れた固有値は、図６に示すように台形積分法でも大きな誤差がのこり、正確に計算するなら、例えばＴａｙｌｏｒ展開法の１０次のような計算アルゴリズム（図７では分かりやすいように８次を使用）を必要とするなど、虚数部において説明したように、そもそもにおいて固有値の絶対値が大きい場合は、前述の（７）式及び（８）式において次数Ｎを相当に大きくしないと収束しないため、台形法のような次数の小さな計算アルゴリズムでは、打切り誤差が大きく、正確な計算ができない。
そのため、“λ＝−３”より小さい固有値を“λ＝−３”に置き換えたことによる差異は、“λ＝−３”のような絶対値の大きな固有値を持つ状態方程式を、台形法のような現在広く知られている計算手法で計算した場合に避け得ない計算誤差を考慮すれば、無視できるものである。更に、“−２”程度を閾値にして、これより小さな固有値をすべて“λ＝−２”に置き換えて、４次のＲｕｎｇｅ−Ｋｕｔｔａ法を使用して“λ＝−２”のときの時間応答を正確に計算したほうが、結果的に妥当な計算結果を得ることができる場合もある。
【０１１５】
なお広く知られたことであるが、固有値が、下記（３３）式で示される任意の複素数である場合を考える。
【数３３】
0006905296.tif 000042
この場合、下記（３４）式が成立する。
【０１１６】
【数３４】
0006905296.tif 000043
従って、複素数の固有値に対する時間応答特性は、図８（ｂ）又は図９（ｂ）に例示した虚数部の時間応答特性と、図１０に例示した実部の時間応答特性との積となり、図１１に例示されるような特性になる。従って、固有値の実部の時間応答特性と虚数部の時間応答特性は、独立に考えることができる。
【０１１７】
次に、本発明の第３の実施形態について説明する。この実施形態は、図３の本発明の第１の実施形態と同様に、ＨＩＬＳ計算装置の実施形態である。図１２は、本発明の第３の実施形態であるＨＩＬＳ計算装置１２００の実施形態を示すブロック図である。図１２において、図３の場合と同じ参照番号を付したブロックは、図３の場合と同じ機能を有する。図１２のＨＩＬＳ計算装置１２００の構成が図３のＨＩＬＳ計算装置３００の構成と異なる点は、まず、図３の固有値実部変更手段３０３及び固有値虚数部変更手段３０４の代わりに、図１２の固有値変更手段１２０１が具備される点である。次に、図１２の行列算出手段１２０２の動作が、図３の行列算出手段３０５の動作とは異なる。更に、図１２の状態方程式計算手段１２０３の動作が図３の状態方程式計算手段３０９の動作と異なる。また、補足的に、図１２のＨＩＬＳ計算装置１２００には、制御装置１２１０が接続される。制御装置１２１０内の特には図示しない第１の制御装置側物理的出力手段及び第２の制御装置側物理的出力手段は、夫々、ＨＩＬＳ計算装置１２００内のデジタル信号入力手段３０６及び電圧入力手段３０７に接続される。デジタル信号入力手段３０６及び電圧入力手段３０７は、計算装置側物理的入力手段として機能する。一方、ＨＩＬＳ計算装置１２００内の電圧出力手段３１１は、制御装置１２１０内の特には図示しない制御装置側物理的入力手段に接続される。電圧出力手段３１１は、計算装置側物理的出力手段として機能する。
【０１１８】
まず、図１２の固有値変更手段１２０１の詳細動作について説明する。前述した図３の本発明の第１の実施形態や図７の本発明の第２の実施形態においては、固有値λの実部σと虚数部ωｉは、夫々、固有値実部変更手段３０３及び固有値虚数部変更手段３０４によって独立に変更可能であった。ここで、Ｒｕｎｇｅ−Ｋｕｔｔａ法やＴａｙｌｏｒ展開法では、図１や図２に示されるように、積λｈの安定領域が概ね半円を描く。従って、これらの数値計算法を使用する場合には、使用する計算アルゴリズムによって、固有値λの実部σ及び虚数部ωｉの間になんらかの拘束条件を設けることが、適切と考えられる。そこで、固有値変更手段１２０１では、使用する数値計算アルゴリズムによって定められる所定の値γ及びμを用いて、実部σ及び虚数部ωの線形結合、即ち下記（３６）式及び（３８）式を満たす関係が設定される。
【０１１９】
即ち、固有値変更手段１２０１は、以下の動作を実行する。まず、虚数部ωｉの係数ωが０以上の場合、実部σ及び虚数部ωｉが、所定の値γ及びμを用いて、下記（３５）式の関係を満たすか否かが判定される。
【数３５】
0006905296.tif 000044
上記（３５）の関係が満たされるならば、下記（３６）式の関係が成立する。
【数３６】
0006905296.tif 000045
上記（３５）式及び（３６）式の関係より、図１２の固有値変更手段１２０１は、固有値λの虚数部ωｉの係数ωが０以上である場合において、（３５）式の関係が満たされる場合には、（３６）式の関係が成立するようにσ又はωの一方又は両方を変更して固有値λを変更する。
一方、虚数部ωｉの係数ωが０未満の場合、実部σ及び虚数部ωｉが、所定の値γ及びμを用いて、下記（３７）式の関係を満たすか否かが判定される。
【数３７】
0006905296.tif 000046
上記（３７）の関係が満たされるならば、下記（３８）式の関係が成立する。
【数３８】
0006905296.tif 000047
上記（３７）式及び（３８）式の関係より、図１２の固有値変更手段１２０１は、固有値λの虚数部ωｉの係数ωが０未満である場合において、（３７）式の関係が満たされる場合には、（３８）式の関係が成立するようにσ又はωの一方又は両方を変更して固有値λを変更する。
【０１２０】
図１３は、第３の実施形態における固有値λの変更方法を説明する図である。図１３に示される複素平面図は、第１の実施形態における図２の図と同じである。図１３の安定領域曲線１３０１、１３０２、１３０３、１３０４、１３０５、及び１３０６は夫々、図２の安定領域曲線２０１、２０２、２０３、２０４、２０５、及び２０６と同じである。
【０１２１】
今、前述した（３６）式は、図１３の複素平面上では、直線１３１０に対応する。また、前述した（３５）式が満たされる領域は、例えば、虚数部の係数ωが図１３の縦軸上で０以上の領域内の上記安定領域曲線１３０１〜１３０６の外側の領域である。例えば座標１３０２が、上記領域に相当する。従って、σとωが、座標１３０２で代表される関係にあるときには、σとωの座標が例えば直線１３１０上の位置１３２１の座標に移動するように、σ又はωの一方の値或いは両方の値が変更される。これにより、変更後の固有値は、安定領域に入ることになる。
【０１２２】
同様に、前述した（３８）式は、図１３の複素平面上では、直線１３１１に対応する。また、前述した（３７）式が満たされる領域は、例えば、虚数部の係数ωが図１３の縦軸上で０未満の領域内の上記安定領域曲線１３０１〜１３０６の外側の領域である。例えば座標１３２２が、上記領域に相当する。従って、σとωが、座標１３２２で代表される関係にあるときには、σとωの座標が例えば直線１３１１上の位置１３２３の座標に移動するように、σ又はωの一方の値或いは両方の値が変更される。これにより、変更後の固有値は、安定領域に入ることになる。
【０１２３】
上記固有値の変更原理に基づく、固有値変更手段１２０１の動作について説明する。図１４は、固有値変更手段１２０１の動作例を示すフローチャートである。
【０１２４】
まず、図３の固有値計算手段３０２から、対角行列Λを取得する（図１４のステップＳ１４０１）。
【０１２５】
次に、ｎ個の固有値λｓ（１≦ｓ≦ｎ）について処理を繰り返すための変数ｓに値１をセットする（図１４のステップＳ１４０２）。
【０１２６】
その後、変数ｓの値が上記ステップＳ１４０２で値１にセットされた後、図１４のステップＳ１４１０で変数ｓの値が＋１ずつインクリメントされながら図１４のステップＳ１４１１で変数ｓの値が値ｎを超えたと判断するまで、図１４のステップＳ１４０３からステップＳ１４０９までの一連の処理を繰り返し実行する。以下、これら一連の処理について順次説明する。
【０１２７】
まず、ステップＳ１４０１で取得した対角行列Λから、変数ｓの値が示す位置の１つの対角要素の固有値“λｓ＝σｓ＋ωｓｉ”を取得する（図１４のステップＳ１４０３）。ここで、σｓは固有値λｓの実部、ωｓｉはｉを虚数単位として固有値λｓの虚数部である。
【０１２８】
次に、ステップＳ１４０３で取得した固有値λｓの虚数部ωｓの値が０以上であるか否かを判定する（図１４のステップＳ１４０４）。
【０１２９】
ステップＳ１４０４の判定がＹＥＳならば、前述した（３５）式に対応する下記（３９）式が満たされるか否かを判定する（図１４のステップＳ１４０５）。
【数３９】
0006905296.tif 000048
【０１３０】
ステップＳ１４０５の判定がＹＥＳならば、前述した（３６）式に対応する下記（４０）式の関係が成立するように、例えば下記（４１）式に従って実部σｓの値を変更する（図１４のステップＳ１４０６）。
【数４０】
0006905296.tif 000049
【数４１】
0006905296.tif 000050
なお、（４０）式の関係が成立するように、実部σｓと虚数部ωｓの両方が変更されてもよい。或いは、（４０）式の関係が成立するように、虚数部ωｓのみが変更されてもよい。
【０１３１】
ステップＳ１４０５の判定がＮＯならば、上記ステップＳ１４０６の変更処理は実行されず、実部σｓと虚数部ωｓは元の値のままとされる。
【０１３２】
前述したステップＳ１４０４の判定がＮＯならば、前述した（３７）式に対応する下記（４２）式が満たされるか否かを判定する（図１４のステップＳ１４０７）。
【数４２】
0006905296.tif 000051
【０１３３】
ステップＳ１４０７の判定がＹＥＳならば、前述した（３８）式に対応する下記（４３）式の関係が成立するように、例えば下記（４４）式に従って実部σｓの値を変更する（図１４のステップＳ１４０８）。
【数４３】
0006905296.tif 000052
【数４４】
0006905296.tif 000053
なお、（４３）式の関係が成立するように、実部σｓと虚数部ωｓの両方が変更されてもよい。或いは、（４３）式の関係が成立するように、虚数部ωｓのみが変更されてもよい。
【０１３４】
ステップＳ１４０７の判定がＮＯならば、上記ステップＳ１４０８の変更処理は実行されず、実部σｓと虚数部ωｓは元の値のままとされる。
【０１３５】
以上のステップＳ１４０６の処理の後、ステップＳ１４０５の判定がＮＯの場合、ステップＳ１４０８の処理の後、又はステップＳ１４０７の判定がＮＯの場合に、適宜変更後の“σｓ＋ωｓｉ”を、修正対角行列（修正固有値包含対角行列）
0006905296.tif 000054
において、変数ｓの現在値が示す位置の対角要素の固有値λｓとして設定する（図１４のステップＳ１４０９）。
【０１３６】
ステップＳ６０８の処理の後、変数ｓの値を＋１インクリメントする（図１４のステップＳ１４１０）。
【０１３７】
そして、インクリメント後の変数ｓの値が対角行列Λの次数ｎを超えたか否かを判定する（図１４のステップＳ１４１１）。
【０１３８】
ステップＳ１４１１の判定がＮＯならば、ステップＳ１４０３の処理に戻って、新たな変数ｓの値が示す対角行列Λの対角要素位置の固有値λｓの変更処理を実行する。
【０１３９】
ステップＳ１４１１の判定がＹＥＳになると、図１２の固有値変更手段１２０１の今回の計算ステップにおける処理を終了し、修正対角行列
0006905296.tif 000055
を、後述する図１２の行列算出手段３０５に出力する。
【０１４０】
以上のようにして、本実施形態では、実行される状態方程式が常に積ｈλに関する安定領域内で解かれるように、対角行列Λの各対角要素の固有値を変更することが可能となる。
【０１４１】
図１２の説明に戻り、次に行列算出手段１２０２について説明する。図３の第１の実施形態では、状態方程式計算手段３０９は、複素固有値を対角要素に持つ対角行列Λを用いて、前述した（３）式及び（４）式に基づいて、解法演算を実行した。しかし、この解法演算は複素演算となるため、対角行列Λの特性によっては、計算時間を要する場合があり得る。
【０１４２】
そこで、第３の実施形態では、行列算出手段１２０２が、固有値変更手段１２０１で算出された修正対角行列
0006905296.tif 000056
から、修正係数行列
0006905296.tif 000057
及び
0006905296.tif 000058
を得る。そして、図１２の状態方程式計算手段１２０３が、これらの修正係数行列に基づいて修正した状態方程式に対して、実数演算による計算を実行できるようにする。
【０１４３】
まず、対角行列Λを得るために使用した前述した（３）式における正則行列Ｐ及びその逆行列Ｐ−１に基づいて、下記（４５）式の変換式が得られる。
【数４５】
0006905296.tif 000059
行列算出手段１２０２は、上記（４５）式の演算を実行する。即ち、図１２において、行列算出手段１２０２は、前述した図１４のフローチャートに基づいて動作する固有値変更手段１２０１から、修正対角行列
0006905296.tif 000060
を入力する。また、行列算出手段１２０２は、固有値計算手段３０２から、前述した正則行列Ｐ及びその逆行列Ｐ−１を入力する。行列算出手段１２０２は、これらの入力したデータに基づいて、（４５）式により、係数行列Ａに対応する修正係数行列
0006905296.tif 000061
を算出し、状態方程式計算手段１２０３に出力する。
なお、行列算出手段１２０２は、状態方程式計算手段１２０３における後述する状態方程式の計算のための補足として、上記修正係数行列の逆行列
0006905296.tif 000062
も算出し、状態方程式計算手段１２０３に出力する。
【０１４４】
また、行列算出手段１２０２は、係数行列Ｂについても修正処理を実行する。図３の第１の実施形態において行列算出手段３０５は、前述した（２４）式の演算により、係数行列Ｂを修正係数行列
0006905296.tif 000063
に修正した。この（２３）式と、前述した（３）式、及び上記（４５）式に基づいて、下記（４６）式の演算式が得られる。
【数４６】
0006905296.tif 000064
行列算出手段１２０２は、上記（４６）式の演算を実行する。即ち、図１２において、行列算出手段１２０２は、通信手段３０１から係数行列Ａを入力し、それに対応する逆行列Ａ−１を算出する。また、行列算出手段１２０２は、通信手段３０１から係数行列Ｂを入力する。行列算出手段１２０２は、これらのデータと、前述の（４５）式により算出した修正係数行列
0006905296.tif 000065
とに基づいて、（４６）式により、修正係数行列
0006905296.tif 000066
を算出し、状態方程式計算手段３０９に出力する。
【０１４５】
最後に、図３の状態方程式計算手段３０９は、下記（４７）式に示される状態方程式を設定する。
【数４７】
0006905296.tif 000067
【０１４６】
状態方程式計算手段３０９は、（４７）式の計算を実際には例えば、前述した（７）式及び（８）式を修正して得られる下記（４８）式及び（４９）式による、例えば５次のＴａｙｌｏｒ展開法の差分方程式として実行する。
【数４８】
0006905296.tif 000068
【数４９】
0006905296.tif 000069
【０１４７】
このようにして、図１２から図１４で説明した第３の実施形態では、固有値変更手段１２０１での固有値の変更を更に適切に実行できると共に、状態方程式計算手段１２０３での状態方程式の解法演算を更に効率的に実行することが可能となる。
【０１４８】
次に、第４の実施形態について説明する。前述した（２）式の計算方法として、例えば前述した（２７）式及び（２８）式においてＮ＝４として４次のＴａｙｌｏｒ展開法を用いた場合を仮定する。この場合、固有値包含対角行列または修正固有値包含対角行列が含む或る固有値λsと或る計算ステップｈに対しては、下記（５０）式及び（５１）式として示される計算が実質的に計算される。
【数５０】
0006905296.tif 000070
【数５１】
0006905296.tif 000071
この時α＝０となる積ｈλｓが存在する。具体的には、ｉを虚数単位として、下記（５２）式又は（５３）式の場合に、上記（５１）式においてα＝０となる。
【数５２】
0006905296.tif 000072
【数５３】
0006905296.tif 000073
の場合に、上記（５１）式において、α＝０となる。この条件を満たすhλｓをゼロ点と言うが、これらのゼロ点を複素平面上で表すと、図１５に示されるごとくとなる。図１５に示される複素平面図は、第１の実施形態における図２の図と同じである。図１５の安定領域曲線１５０１、１５０２、１５０３、及び１５０４は夫々、図２の安定領域曲線２０１、２０２、２０３、及び２０４と同じである。図１５において、座標点１５１０及び１５１１が（５２）式の場合に対応し、座標点１５１２及び１５１３が（５３）式の場合に対応する。図１５から理解されるように、これらのゼロ点は、４次のＴａｙｌｏｒ展開法の安定領域曲線１５００３の内側、即ち安定領域内に存在することがわかる。
【０１４９】
このような固有値の場合、前述した（５０）式は、下記（５４）式のようになる。
【数５４】
0006905296.tif 000074
この値はｕｎが一定であれば、ｘ∞の値、つまり最終値であり、最初の１ステップで最終値に収束してしまい、過渡的な値を得ることができない。最終結果を得るだけの計算であればよいが、ＳＩＬＳやＨＩＬＳの場合過渡応答も極めて重要であるため、この状態は好ましいものではない。
【０１５０】
そこで、第４の実施形態では、このような固有値は、例え計算アルゴリズムの安定領域内であっても、他の値に変更する。第４の実施形態では、前述した図１２の第３の実施形態の場合と同様の構成において、固有値変更手段１２０１が、（５１）式においてα＝０となる条件を充足することを判別するゼロ点検出手段を更に備える。そして、このゼロ点検出手段は、上記条件が充足されたことを判別した場合に、固有値λｓを他の固有値に変更する、例えば、下記（５５）式のように固有値の変更を行う。あるいは、５次のＴａｙｌｏｒ展開法で計算する。
【数５５】
0006905296.tif 000075
【０１５１】
前述した（２）式において、下記（５６）式の場合には、その固有値は前述した（５２）式の値となり、前述した（５１）式に対する判別はα＝０となる。
【数５６】
0006905296.tif 000076
【０１５２】
図１６は、係数行列Ａ、Ｂを持つ前述した（２）式を、下記の３つの条件で計算した場合における、状態ベクトルＸの第１要素ｘ１と第２要素ｘ２の時系列計算結果を示す図である。
１６０１：オリジナルの４次のＴａｙｌｏｒ展開法で計算した場合
１６０２：前述の（５５）式のように固有値の変更を行った上で、
４次のＴａｙｌｏｒ展開法で計算した場合
１６０３：５次のＴａｙｌｏｒ展開法で計算した場合
また、図１６には、理論値のオリジナルの時系列計算結果１６０４と、理論値×０．７５の時系列計算結果１６０５も、一緒に示されている。
図１６では、２次方程式
【数５７】
0006905296.tif 000077
の固有値が、図１２の第３の実施形態の場合と同様の固有値変更手段１２０１により、前述の（５２）式の場合の値とされる。
次に、図１２の第３の実施形態の場合と同様の行列算出手段１２０２により、修正係数行列
0006905296.tif 000078
及び
0006905296.tif 000079
が生成される。
そして、図１２の第３の実施形態の場合と同様の状態方程式計算手段１２０３で前述の（４７）式の状態方程式が計算される。図１６（ａ）は横軸：時間（秒）に対する縦軸：ｄｘ／ｄｔの計算結果を示しており、図１６（ｂ）は横軸：時間（秒）に対する縦軸：ｘの計算結果を示している。
【０１５３】
図１６において、固有値が前述した（５２）式として示される値のままの計算結果１６０１では、振動現象は表現できていない。これに対して、第４の実施形態において固有値が０．７５倍された計算結果１６０２では、真値１６０４とは異なるものの、振動現象は良く表現されていることがわかる。５次のＴａｙｌｏｒ展開法の計算結果１６０３も、やはり真値１６０４とは異なる。
計算結果が真値と異なるのは、固有値の絶対値が大きく４次や５次では計算誤差が大きいためであるので、例えば１０次のＴａｙｌｏｒ展開法などを使用すれば、真値に近い値を得ることができる。
固有値が前述した（５３）式として示される値も、第４の実施形態は同様に適用できる。
【０１５４】
上述の第４の実施形態では、固有値変更手段１２０１内のゼロ点検出手段が条件が充足されたことを判別した場合に、例えば下記（５５）式のように、固有値λｓを他の固有値に変更した。これに対し、ゼロ点検出手段が条件が充足されたことを判別した場合に、固有値を変更するのではなく、計算次数Ｎを変更するようにしてもよい。
【０１５５】
以上説明したように、固有値計算手段３０２と、固有値実部変更手段３０３／固有値虚数部変更手段３０４又は固有値変更手段１２０１とを有する上述の各実施形態によるＨＩＬＳ，ＳＩＬＳ計算装置は、いわゆるＳｔｉｆｆ状態方程式であっても、計算精度が高く、計算速度の速い数値計算手法を用いることが可能となる。
これらの構成がＳＩＬＳに適用された場合は、計算精度が向上すると同時に計算時間を短縮することが可能となる。
上記構成がＨＩＬＳに適用された場合には、従来では不可能であったＳｔｉｆｆ状態方程式で表される制御対象を、リアルタイム処理として模擬することが可能となる。
この結果、制御装置の開発を飛躍的に効率化することが可能となり、産業界の発展に寄与するものである。
【符号の説明】
【０１５６】
１０１、１０２、１０３、２０１、２０２、２０３，２０４、２０５、２０６、１３０１、１３０２、１３０３、１３０４、１３０５、１３０６、１５０１、１５０２、１５０３、１５０４、  安定領域曲線
２０７  ナイキスト周波数を超えない範囲
３００、１２００  ＨＩＬＳ計算装置
３０１  通信手段
３０２  固有値計算手段
３０３  固有値実部変更手段
３０４  固有値虚数部変更手段
３０５  行列算出手段
３０６  デジタル信号入力手段
３０７  電圧入力手段
３０８  ベクトル算出手段
３０９  状態方程式計算手段
３１０  時計装置
３１１  電圧出力手段
３２０  ホストコンピュータ
７００  ＳＩＬＳ計算装置
７０１  制御装置制御則計算手段
７０２  行列設定手段
７０３  ディスプレイ装置
【要約】
【課題】制御対象の動作を数学的にモデル化した状態方程式を解く演算を実行することにより制御対象の動作を模擬し、それにより制御装置の設計及び検証を支援する計算装置の技術に関し、状態方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とする。
【解決手段】固有値計算手段３０２は、変数ベクトル項の係数行列を、対角要素が固有値である対角行列又はジョルダン標準形行列の何れかの行列である固有値包含対角行列に変換する。固有値変更手段１２０１は、固有値包含対角行列内の各固有値が所定の条件を満たす場合に、固有値の実部又は虚数部の一方又は両方を所定の定数に変更し、修正固有値包含対角行列を出力する。方程式変換手段である行列算出手段１２０２は、線形状態方程式を修正固有値包含対角行列に基づく形式に変換し、修正線形状態方程式を出力する。状態方程式計算手段１２０３は、修正線形状態方程式を差分方程式の形式に変換し、差分方程式に対する解法演算を実行することにより、制御装置への出力信号を生成する。
【選択図】図１２【図１】
0006905296.tif 000128
【図２】
0006905296.tif 000129
【図３】
0006905296.tif 000130
【図４】
0006905296.tif 000131
【図５】
0006905296.tif 000132
【図６】
0006905296.tif 000133
【図７】
0006905296.tif 000134
【図８】
0006905296.tif 000135
【図９】
0006905296.tif 000136
【図１０】
0006905296.tif 000137
【図１１】
0006905296.tif 000138
【図１２】
0006905296.tif 000139
【図１３】
0006905296.tif 000140
【図１４】
0006905296.tif 000141
【図１５】
0006905296.tif 000142
【図１６】
0006905296.tif 000143