(57)【特許請求の範囲】
【請求項１】
計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する回路方程式作成部と、
前記回路方程式作成部で作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う順序付け部と、
前記順序付け部で直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行うＬＵ分解部と、
前記ＬＵ分解部でＬＵ分解された行列を用いて求解演算を行う求解演算部と、
を備える周波数特性解析装置。
【請求項２】
前記順序付け部は、前記回路方程式作成部で最初に回路方程式が作成された場合に、前記完全ピボッティングによる順序付けを行う、請求項１記載の周波数特性解析装置。
【請求項３】
前記ＬＵ分解部でＬＵ分解された行列を用いて計算精度を評価する計算精度評価部を更に備えており、
前記順序付け部は、前記計算精度評価部で計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う、請求項１又は請求項２記載の周波数特性解析装置。
【請求項４】
前記計算精度評価部は、前記ＬＵ分解部でＬＵ分解された行列を用いて計算精度の指標を示す不安定因子を求め、該不安定因子が予め規定された閾値を超えた場合に、計算精度の悪化が生じたと評価する、請求項３記載の周波数特性解析装置。
【請求項５】
回路方程式作成部が、計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する第１ステップと、
順序付け部が、前記回路方程式作成部で作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う第２ステップと、
ＬＵ分解部が、前記順序付け部で直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行う第３ステップと、
求解演算部が、前記ＬＵ分解部でＬＵ分解された行列を用いて求解演算を行う第４ステップと、
演算処理装置が、前記計算周波数を変えて前記第１ステップから第４ステップを繰り返す第５ステップと、
を有する周波数特性解析方法。
【請求項６】
計算精度評価部が、前記ＬＵ分解部でＬＵ分解された行列を用いて計算精度を評価する第６ステップを更に有しており、
前記順序付け部は、前記回路方程式作成部で最初に回路方程式が作成された場合、或いは前記計算精度評価部で計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う、請求項５記載の周波数特性解析方法。
【請求項７】
コンピュータを、
計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する回路方程式作成手段と、
前記回路方程式作成手段で作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う順序付け手段と、
前記順序付け手段で直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行うＬＵ分解手段と、
前記ＬＵ分解手段でＬＵ分解された行列を用いて求解演算を行う求解演算手段と、
して機能させる周波数特性解析プログラム。
【請求項８】
前記コンピュータを、前記ＬＵ分解手段でＬＵ分解された行列を用いて計算精度を評価する計算精度評価手段として更に機能させ、
前記順序付け手段は、前記回路方程式作成手段で最初に回路方程式が作成された場合、或いは前記計算精度評価手段で計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う、請求項７記載の周波数特性解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、周波数特性解析装置、方法、及びプログラムに係り、特に電力系統の周波数特性解析等に適用して有用な周波数特性解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
電力系統で発生する過電圧や過電流の解析、電力品質の解析、交直変換器の動特性の解析等には、電力及び電流の時間変化を波形レベルで解析可能な瞬時値解析が有効である。このような瞬時値解析が可能な瞬時値解析プログラムの多くには、特定の母線から系統を見込んだインピーダンスの周波数特性を計算する機能があり、直流送電設備や周波数変換設備の交直変換器と交流系統の相互作用によって発生する高調波不安定現象や高調波特性を含めた過渡過電圧の解析等で用いられる。
【０００３】
図１２は、従来の周波数特性計算の計算過程を示すフローチャートである。図１２に示す通り、周波数特性計算は、通常、回路の交流定常解を求める計算（交流定常解計算）を、周波数を変化させながら繰り返し実行する計算である。交流定常解計算は、連立一次方程式となる回路方程式を作成する定式化過程と、作成された回路方程式を解く求解過程とからなる。
【０００４】
定式化過程は、回路中の全ノードと全回路部品の情報から回路方程式を作成する過程である。例えば、回路にＮｎ個の独立ノード、Ｎｂ個の回路部品（ブランチ）が存在すると仮定した場合、回路中のＮｎ個のノードの電圧からなるベクトルｕ、Ｎｂ個の回路部品を流れる電流からなるベクトルｉ、Ｎｂ個の回路部品の両端電圧からなるベクトルｖを未知変数として取り扱い、その回路方程式は以下の（１）式となる。
【数１】
0006650333.tif 000002
【０００５】
上記（１）式をスパースタブロー方程式と呼び、Ｆをタブロー行列と呼ぶ。上記（１）式中のＦ，ｘ，ｙは、以下の（２）式で表される。
【数２】
0006650333.tif 000003
【０００６】
ここで、Ａは、回路中のＮｎ個の独立ノードとＮｂ個のブランチの接続状態を関連付けるＮｂ×Ｎｎのブランチ対ノード接続行列である。Ｉは、Ｎｂ×Ｎｂの単位行列である。Ｂ＝［Ｂ１ −Ｂ２］は、各ブランチのブランチ電圧ｖとブランチ電流ｉとを関連付けるＮｂ×２Ｎｂのブランチ特性行列である。ｓは、回路中の電圧源、電流源の値により定まるベクトルである。Ｎ＝Ｎｎ×２Ｎｂのとき、ｘ，ｙは要素数がＮ個のベクトルとなり、ＦはＮ×Ｎの行列となる。交流定常解計算では、ここで示したスパースタブロー方程式の全ての量は複素数となる。
【０００７】
ブランチ特性行列Ｂは、ブランチ特性方程式から得られる。周波数をｆ＝ω／（２π）としたとき、例えばインダクタＬ、キャパシタＣのブランチ特性方程式はそれぞれ、以下の（３）式及び（４）式で表される。
【０００８】
【数３】
0006650333.tif 000004
【数４】
0006650333.tif 000005
【０００９】
但し、ｉＬ，ｉＣは、Ｌ，Ｃに流れる電流（ブランチ電流）であり、ｖＬ，ｖＣは、Ｌ，Ｃの両端電圧（ブランチ電圧）である。Ｌがｂ番目のブランチだった場合、Ｂ及びｓのｂ行目は、以下の（５）式で表される。
【数５】
0006650333.tif 000006
【００１０】
但し、上記（５）式左辺の「｜」印より左は行列Ｂ１のｂ行目を示し、「｜」印より右は行列−Ｂ２のｂ行目を示し、右辺はｓのｂ行目を示す。また、上記（５）式下部の２つの矢印は、１と−１／ｊωＬがそれぞれＢ１とＢ２のｂ列目であることを示し、「…」で省略された要素は０である。尚、解析回路中に非線形素子が配置された場合には、非線形素子を無視するか、或いは原点近くの特性で線形化して取り扱う。また、回路中の電圧源、電流源の値により定まるベクトルｓは、電圧源が発生する電圧波形がＥｍｃｏｓ（ωｔ＋θ）のとき、ｓの該当する要素にＥｍｅｊθを代入する。Ｉｍｃｏｓ（ωｔ＋θ）を発生する電流源については、Ｉｍｅｊθを代入する。
【００１１】
求解過程は、順序付け、ＬＵ分解、前進・後退代入と呼ばれる処理が順番に実行される。順序付けは、計算精度の向上を目的に行列要素を並べ替える操作である。例えば、以下に示す手順で行列要素の並べ替えが行われる。
（イ）行列Ｆについて、各行に存在する非零要素数ＮＲと、各列に存在する非零要素数ＮＣとを求める。
（ロ）ＮＲＣ＝（ＮＲ−１）×（ＮＣ−１）を要素毎に計算し、ＮＲＣが最も小さい要素を１行目の対角要素（ピボット）となるように入れ替え操作を行う。ＮＲＣが同じ要素が複数ある場合は、最も絶対値が大きい要素を選択する。
（ハ）１行１列を除いた残りの行列（部分行列）をＦとして考え、（イ）〜（ハ）の手順を最後の行まで繰り返し行う。
【００１２】
ここで、手順（ロ）は、ＮＲＣが最も小さい要素（ＮＲＣが同じ場合は絶対値が大きい要素）を選び出し、それを対角要素にする作業であるが、例えばｉ行目の対角要素の探索において、その探索対象をｉ行の要素のみとした場合と、既に対角要素を決定した行と列を除く全ての要素とした場合の２通りの手法が考えられる。前者は部分ピボッティング、後者は完全ピボッティングと呼ばれている。一般に、部分ピボッティングは、探索範囲が狭いため計算量は少なく、完全ピボッティングは、探索範囲が広いため計算量は大きいが、より最適な対角要素を選び出すことができるため計算精度は高くなる。
【００１３】
ＬＵ分解及び前進・後退代入は、回路方程式の係数行列（上記（１）式中のタブロー行列Ｆ）を、下三角行列Ｌと上三角行列Ｕとに分解し、前進代入及び後退代入と呼ばれる操作によって未知数となるｘを求める操作である。いま、上記（１）式の行列要素をａｉｊ（但し、１番目の添え字は行番号、２番目の添え字は列番号、ｉ＝１，２，…，Ｎ、ｊ＝１，２，…，Ｎ）としたとき、ＦのＬＵ分解後の下三角行列Ｌの要素ｌｉｊと上三角行列Ｕの要素ｕｉｊはそれぞれ以下の（６）式及び（７）式により計算することができる。
【００１４】
【数６】
0006650333.tif 000007
【００１５】
【数７】
0006650333.tif 000008
【００１６】
次に、求めたｕｉｊ，ｌｉｊを用いた前進・後退代入によりｘを求める。尚、ＬＵ分解に続く前進・後退代入の過程の詳細な説明は省略する。ｘを得ることにより、交流定常解として前述したノード電圧、ブランチ電流、ブランチ電圧が求められる。以上説明した一連の過程（定式化、求解）を、計算周波数ｆを更新して繰り返し行うことで、回路の周波数特性を求めることができる。尚、以下の非特許文献１，２には、周波数特性計算に係るものではないが、電圧及び電流の時間変化を求める通常の瞬時値解析における計算時間の短縮を図る方法が開示されている。
【先行技術文献】
【非特許文献】
【００１７】
【非特許文献１】米澤力道，野田琢，「電力系統瞬時値解析プログラムの開発（その５）−大規模疎行列ＬＵ分解アルゴリズムの改良によるＸＴＡＰの高速化−」，電力中央研究所報告  研究報告Ｈ１１００９，２０１２
【非特許文献２】米澤力道，野田琢，「電力系統瞬時値解析プログラムの開発（その７）−回路部品の内部状態更新処理の並列化によるＸＴＡＰの高速化−」、電力中央研究所報告  研究報告Ｈ１３００５，２０１４
【発明の概要】
【発明が解決しようとする課題】
【００１８】
ところで、上述したＬＵ分解に基づく求解は、周波数特性計算で行われるだけでなく、電圧及び電流の時間変化を求める通常の瞬時値解析でも行われる。通常の瞬時値解析では、上記（１）式を各計算時間ステップで解く計算過程となるが、素子の値が時間によって変化する時変素子や非線形素子が回路中に存在しない場合には、上記（１）式で計算時間ステップごとに変化するのはｙのみであり、Ｆは変化しない。即ち、通常の瞬時値解析では、Ｆに対する順序付け及びＬＵ分解を最初に一度だけ行っておけば良く、各計算時間ステップの解ｘは、変化するｙに対して前進代入及び後退代入を行うだけで求めることができる。
【００１９】
しかしながら、周波数特性計算においては、計算周波数ｆの更新によってブランチ特性方程式が変化するため、設定された計算周波数ｆ毎にＦが変化してしまう。このため、周波数特性計算における求解過程では、通常の瞬時値解析における求解過程のように順序付け及びＬＵ分解を省略することはできず、計算周波数ｆ毎に順序付け及びＬＵ分解を実行する必要がある。このように、周波数特性計算では、瞬時値解析の場合とは異なり、求解過程により多くの時間を要し、その計算時間の短縮が課題となっている。
【００２０】
ここで、前述の通り、求解過程における順序付けでは、部分ピボッティング又は完全ピボッティングを用いることができる。周波数特性計算における求解過程の計算時間を短縮するには、完全ピボッティングよりも計算時間が少ない部分ピボッティングを用いれば良いとも考えられる。部分ピボッティングは、計算精度は原理的に完全ピボッティングより劣るものの、多くの場合で解析精度の低下は問題とならないが、部分ピボッティングで計算精度に問題を生じる例が報告されており、必ずしも万能の手法とはいえない。そこで、計算精度の維持の観点から完全ピボッティングを用いつつ、求解過程の計算時間を短縮する必要がある。
【００２１】
本発明は上記事情に鑑みてなされたものであり、計算精度を維持しつつ周波数特性計算に要する時間を大幅に短縮することが可能な周波数特性解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００２２】
上記課題を解決するために、本発明の周波数特性解析装置（１）は、計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する回路方程式作成部（１４ａ）と、前記回路方程式作成部で作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う順序付け部（１４ｂ）と、前記順序付け部で直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行うＬＵ分解部（１４ｃ）と、前記ＬＵ分解部でＬＵ分解された行列を用いて求解演算を行う求解演算部（１４ｅ）と、を備える。
また、本発明の周波数特性解析装置は、前記順序付け部が、前記回路方程式作成部で最初に回路方程式が作成された場合に、前記完全ピボッティングによる順序付けを行う。
また、本発明の周波数特性解析装置は、前記ＬＵ分解部でＬＵ分解された行列を用いて計算精度を評価する計算精度評価部（１４ｄ）を更に備えており、前記順序付け部が、計算精度評価部で計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う。
また、本発明の周波数特性解析装置は、前記計算精度評価部が、前記ＬＵ分解部でＬＵ分解された行列を用いて計算精度の指標を示す不安定因子（ρ）を求め、該不安定因子が予め規定された閾値（ρｔｈ）を超えた場合に、計算精度の悪化が生じたと評価する。
本発明の周波数特性解析方法は、計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する第１ステップ（Ｓ３０）と、前記第１ステップで作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う第２ステップ（Ｓ４１）と、前記第２ステップで直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行う第３ステップ（Ｓ４２）と、前記第３ステップでＬＵ分解された行列を用いて求解演算を行う第４ステップ（Ｓ４４）と、前記計算周波数を変えて前記第１ステップから第４ステップを繰り返す第５ステップ（Ｓ２０）と、を有する。
また、本発明の周波数特性解析方法は、前記第３ステップでＬＵ分解された行列を用いて計算精度を評価する第６ステップ（Ｓ４３）を更に有しており、前記第２ステップが、前記第１ステップで最初に回路方程式が作成された場合、或いは前記第６ステップで計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う。
本発明の周波数特性解析プログラムは、コンピュータを、計算周波数に応じた解析対象の回路方程式を、疎行列を用いて作成する回路方程式作成手段（１４ａ）と、前記回路方程式作成手段で作成される回路方程式のうち、特定の場合に作成された回路方程式の前記疎行列に対してのみ完全ピボッティングによる順序付けを行う順序付け手段（１４ｂ）と、前記順序付け手段で直近に行われた順序付け結果を用いて前記疎行列のＬＵ分解を行うＬＵ分解手段（１４ｃ）と、前記ＬＵ分解手段でＬＵ分解された行列を用いて求解演算を行う求解演算手段（１４ｅ）と、して機能させる。
また、本発明の周波数特性解析プログラムは、前記コンピュータを、前記ＬＵ分解手段でＬＵ分解された行列を用いて計算精度を評価する計算精度評価手段（１４ｄ）として更に機能させ、前記順序付け手段が、前記回路方程式作成手段で最初に回路方程式が作成された場合、或いは前記計算精度評価手段で計算精度の悪化が生じたと評価された場合に、前記完全ピボッティングによる順序付けを行う。
【発明の効果】
【００２３】
本発明によれば、疎行列を用いて計算周波数に応じた解析対象の回路方程式を作成し、作成された回路方程式のうち、特定の場合に作成された回路方程式の疎行列に対してのみ完全ピボッティングによる順序付けを行い、直近に行われた順序付け結果を用いて疎行列のＬＵ分解を行い、ＬＵ分解された行列を用いて求解演算を行うようにしている。これにより、特定の場合以外における完全ピボッティングによる順序付けを省略することができることから、計算精度を維持しつつ周波数特性計算に要する時間を大幅に短縮することが可能であるという効果がある。
【図面の簡単な説明】
【００２４】
【図１】本発明の一実施形態による周波数特性解析装置の要部構成を示すブロック図である。
【図２】本発明の一実施形態による周波数特性解析方法の一例を示すフローチャートである。
【図３】回路規模と計算速度との関係を検証するために用いた検証用回路を示す図である。
【図４】図３中の発電機Ｇ２〜Ｇ１０の模擬に用いた等価回路を示す図である。
【図５】図３中の変圧器ＴＲ２〜ＴＲ１０の模擬に用いた等価回路を示す図である。
【図６】図３に示す「ケース１」〜「ケース９」について周波数特性計算を行った場合に要した計算時間を示す図である。
【図７】回路規模と速度向上倍数との関係を示す図である。
【図８】図３に示す「ケース９」について周波数特性計算を行った場合の各サンプル点における不安定因子を示す図である。
【図９】計算精度及び計算速度を検証するために用いた検証用回路を示す図である。
【図１０】図９の検証用回路に対して行った周波数特性計算の計算結果を示す図である。
【図１１】図９の検証用回路について周波数特性計算を行った場合の各サンプル点における不安定因子を示す図である。
【図１２】従来の周波数特性計算の計算過程を示すフローチャートである。
【発明を実施するための形態】
【００２５】
以下、図面を参照して本発明の一実施形態による周波数特性解析装置、方法、及びプログラムについて詳細に説明する。尚、以下では、解析対象が電力系統である場合を例に挙げて説明するが、本発明は、解析対象が電力系統に限られる訳ではなく、任意の電気回路を解析対象とすることができる。
【００２６】
〈周波数特性解析装置〉
図１は、本発明の一実施形態による周波数特性解析装置の要部構成を示すブロック図である。図１に示す通り、本実施形態の周波数特性解析装置１は、入力装置１１、出力装置１２、補助記憶装置１３、及び演算処理装置１４を備えており、解析対象としての電力系統の周波数特性計算を行う。このような周波数特性解析装置１は、例えばデスクトップ型のパーソナルコンピュータによって実現される。
【００２７】
入力装置１１は、例えばキーボードやポインティングデバイス等を備えており、周波数特性解析装置１を使用するユーザの操作に応じた指示を演算処理装置１４に出力する。出力装置１２は、例えば液晶表示装置やプリンタ等を備えており、演算処理装置１４から出力される各種情報を出力する。補助記憶装置１３は、例えばＨＤＤ（ハードディスクドライブ）やＳＳＤ（ソリッドステートドライブ）等を備えており、電力系統の解析を行う上で用いられる各種情報（順序付け結果ＯＲ及び評価パラメータＰＲ：詳細は後述する）及び周波数特性解析プログラムＰＧ等を記憶する。
【００２８】
演算処理装置１４は、入力装置１１から入力される操作指示に基づき、補助記憶装置１３に記憶された周波数特性解析プログラムＰＧに従い補助記憶装置１３に記憶された各種情報を用いて電力系統の周波数特性計算を行う。この演算処理装置１４は、回路方程式作成部１４ａ（回路方程式作成手段）、順序付け部１４ｂ（順序付け手段）、ＬＵ分解部１４ｃ（ＬＵ分解手段）、計算精度評価部１４ｄ（計算精度評価手段）、及び求解演算部１４ｅ（求解演算手段）を備える。
【００２９】
回路方程式作成部１４ａは、解析対象である電力系統の回路方程式を、疎行列を用いて作成する。具体的に、回路方程式作成部１４ａは、前述した（１）式に示されるスパースタブロー方程式を、計算周波数ｆに応じて作成する。尚、上記の疎行列は、前述した（２）式に示すタブロー行列Ｆである。解析対象である電力系統の構成が変わらない場合において、計算周波数ｆが変わると、前述した（２）式に示すタブロー行列Ｆ中のブランチ特性行列Ｂ＝［Ｂ１ −Ｂ２］が変わる。
【００３０】
順序付け部１４ｂは、回路方程式作成部１４ａで作成される回路方程式のタブロー行列Ｆに対して、完全ピボッティングによる順序付けを行う。ここで、順序付け部１４ｂは、回路方程式作成部１４ａで作成される回路方程式のうち、特定の場合に作成された回路方程式のタブロー行列Ｆに対してのみ完全ピボッティングによる順序付けを行う。これは、周波数特性解析に要する時間を短縮するためである。
【００３１】
具体的に、順序付け部１４ｂは、回路方程式作成部１４ａで最初に回路方程式が作成された場合、或いは計算精度評価部１４ｄで計算精度の悪化が生じたと評価された場合に、完全ピボッティングによる順序付けを行う。尚、順序付け部１４ｂにより行われた完全ピボッティングによる順序付けの結果は、順序付け結果ＯＲとして補助記憶装置１３に記憶される。
【００３２】
ＬＵ分解部１４ｃは、補助記憶装置１３に記憶された順序付け結果ＯＲを用いてタブロー行列ＦのＬＵ分解を行う。ここで、補助記憶装置１３に記憶された順序付け結果ＯＲは、順序付け部１４ｂによって完全ピボッティングによる順序付けが行われる度に更新される。このため、ＬＵ分解部１４ｃは、順序付け部１４ｂで直近に行われた順序付け結果を用いてタブロー行列ＦのＬＵ分解を行っていると言うことができる。
【００３３】
計算精度評価部１４ｄは、ＬＵ分解部１４ｃでＬＵ分解された行列を用いて計算精度を評価する。具体的に、計算精度評価部１４ｄは、以下の（８）式に示す計算精度の指標を示す不安定因子ρを求め、この不安定因子ρが予め規定された閾値ρｔｈ以上になった場合に、計算精度の悪化が生じたと評価する。尚、以下の（８）式では、下三角行列Ｌの要素ｌｉｊと上三角行列Ｕの要素ｕｉｊをまとめてｌｕｉｊと表記している。また、上記の閾値ρｔｈは、評価パラメータＰＲとして補助記憶装置１３に記憶されている。
【００３４】
【数８】
0006650333.tif 000009
【００３５】
上記の（８）式から明らかな通り、不安定因子ρは、ＬＵ分解前後における行列要素の無限大ノルムの増加比となる。通常、スケーリング（或いは、正規化、均等化）と呼ばれる前処理によって、上記（８）式の分母は１付近の値となるため、不安定因子ρは、概ねlｕｉｊの最大要素の値と等しくなる。順序付けにおいて、絶対値が小さい値等、対角要素として適切ではない要素が対角要素に選択されている場合には、前述した（７）式の分母が小さな値となるため、ｌｉｊ，ｕｉｊが増大し、不安定因子ρが増大する。そのため、不安定因子ρを用いることで順序付けが適切かどうかを判断することが可能である。
【００３６】
求解演算部１４ｅは、ＬＵ分解部１４ｃでＬＵ分解された行列（下三角行列Ｌ、上三角行列Ｕ）を用いて求解演算を行う。具体的には、上三角行列Ｕの要素ｕｉｊ及び下三角行列Ｌの要素ｌｉｊを用いた前進・後退代入によって、前述した（１），（２）式に示すｘを求める。
【００３７】
周波数特性解析プログラムＰＧは、例えばＣＤ−ＲＯＭ又はＤＶＤ（登録商標）−ＲＯＭ等のコンピュータ読み取り可能な記録媒体に記録されたもの、或いはインターネット等のネットワークを介してダウンロードされたものが周波数特性解析装置１にインストールされることにより、補助記憶装置１３に記憶される。補助記憶装置１３に記憶された周波数特性解析プログラムＰＧが読み出されて実行されることにより、周波数特性解析装置１の各ブロックの機能（例えば、回路方程式作成部１４ａ、順序付け部１４ｂ、ＬＵ分解部１４ｃ、計算精度評価部１４ｄ、及び求解演算部１４ｅ）がソフトウェア的に実現される。つまり、これらの機能は、ソフトウェアとハードウェア資源とが協働することによって実現される。
【００３８】
〈周波数特性解析方法〉
次に、上述した周波数特性解析装置１を用いた周波数特性解析方法について説明する。周波数特性の解析を行う場合には、まず、ユーザによって解析対象である電力系統の入力が行われる。具体的には、ユーザによって周波数特性解析装置１の入力装置１１が操作され、ユーザの指示に応じて、電力系統を構成する全てのノード及び回路部品の情報を演算処理装置１４に入力する処理が行われる。
【００３９】
次に、ユーザによって解析条件の入力が行われる。例えば、周波数特性解析装置１の入力装置１１を操作するユーザによって、計算開始周波数ｆｂｅｇｉｎ、計算終了周波数ｆｅｎｄ、サンプル点数Ｋ、及び計算周波数の刻み方（等間隔、或いは対数軸上等間隔）等を入力する処理が行われる。
【００４０】
以上の入力が完了し、ユーザによって解析開始の指示がなされると、図２に示すフローチャートに従って、電力系統の周波数特性計算が開始される。図２は、本発明の一実施形態による周波数特性解析方法の一例を示すフローチャートである。周波数特性計算が開始されると、まず図１に示す演算処理装置１４によって、各種の初期処理行われる（ステップＳ１０）。例えば、サンプル点を特定する変数ｋを初期化し（値を１に設定し）、各サンプル点における計算周波数ｆｋを求める処理が行われる。
【００４１】
次に、演算処理装置１４によって、計算周波数ｆを更新する処理が行われる（ステップＳ２０）。ここでは、上記のステップＳ１０の処理で変数ｋの値が１に設定されていることから、計算周波数ｆは計算周波数ｆ１に更新（設定）される。尚、この計算周波数ｆ１は、例えばユーザによって入力された計算開始周波数ｆｂｅｇｉｎである。これにより、最初のサンプル点が特定される。
【００４２】
次いで、定式化過程（電力系統の回路方程式の作成）が行われる（ステップＳ３０：第１ステップ）。具体的には、ユーザによって入力された電力系統を構成する全てのノード及び回路部品の情報を用い、ステップＳ２０で更新（設定）された計算周波数ｆに応じたスパースタブロー方程式（前述した（１），（２）式に示される方程式）を作成する処理が回路方程式作成部１４ａで行われる。定式化過程が完了すると、続いて求解過程が行われる（ステップＳ４０）。
【００４３】
求解過程が開始されると、まず、回路方程式作成部１４ａで作成された回路方程式のタブロー行列Ｆに対し、完全ピボッティングによる順序付けが順序付け部１４ｂで行われる（ステップＳ４１：第２ステップ）。尚、順序付け部１４ｂで行われた完全ピボッティングによる順序付けの結果は、順序付け結果ＯＲとして補助記憶装置１３に記憶される。次に、補助記憶装置１３に記憶された順序付け結果ＯＲを用いてタブロー行列ＦのＬＵ分解を行う処理がＬＵ分解部１４ｃで行われる（ステップＳ４２：第３ステップ）。つまり、ＬＵ分解部１４ｃの初回の処理（計算）は、順序付け部１４ｂによる順序付けが実施されてから行われる。
【００４４】
次いで、ＬＵ分解部１４ｃでＬＵ分解された行列を用いて計算精度を評価する処理が計算精度評価部１４ｄで行われる（ステップＳ４３：第６ステップ）。具体的には、ＬＵ分解部１４ｃでＬＵ分解された行列を用いて、前述した（８）式に示す計算精度の指標を示す不安定因子ρを求め、この不安定因子ρが閾値ρｔｈよりも小であるか否かを判断する処理が計算精度評価部１４ｄで行われる。尚、上記の閾値ρｔｈは、評価パラメータＰＲとして補助記憶装置１３に記憶されているものである。
【００４５】
不安定因子ρが閾値ρｔｈよりも小であると判断された場合（ステップＳ４３の判断結果が「ｙｅｓ」である場合）には、ＬＵ分解部１４ｃでＬＵ分解された行列を用いた求解演算が求解演算部１４ｅで行われる（ステップＳ４４：第４ステップ）。具体的には、上三角行列Ｕの要素ｕｉｊ及び下三角行列Ｌの要素ｌｉｊを用いた前進・後退代入によって、前述した（１），（２）式に示すｘを求める演算が求解演算部１４ｅで行われる。
【００４６】
以上の処理が終了すると、変数ｋの値がサンプル点数Ｋよりも大になったか否かが、演算処理装置１４によって判断される（ステップＳ５０）。変数ｋの値がサンプル点数Ｋ以下であると判断された場合（ステップＳ５０の判断結果が「ｎｏ」の場合）には、変数ｋの値をインクリメントしてから計算周波数ｆを更新する処理が演算処理装置１４によって行われる（ステップＳ２０：第５ステップ）。ここでは、計算周波数ｆが、計算周波数ｆ２に更新される。これにより、次のサンプル点が特定される。
【００４７】
計算周波数ｆが更新されると、その更新された計算周波数ｆに応じた回路方程式を作成する処理が回路方程式作成部１４ａで行われる（ステップＳ３０）。その後、順序付け部１４ｂによる順序付けが行われることなく、補助記憶装置１３に記憶された順序付け結果ＯＲを用いてタブロー行列ＦのＬＵ分解を行う処理がＬＵ分解部１４ｃで行われる（ステップＳ４２）。つまり、ＬＵ分解部１４ｃで行われる処理（計算）は２回目以降の計算に該当することから、順序付け部１４ｂで直近に行われた順序付け結果を用いてタブロー行列ＦのＬＵ分解が行われる。
【００４８】
続いて、計算精度を評価する処理が計算精度評価部１４ｄで行われる（ステップＳ４３）。そして、不安定因子ρが閾値ρｔｈよりも小であると判断された場合（ステップＳ４３の判断結果が「ｙｅｓ」である場合）には、求解演算（ステップＳ４４）、及び変数ｋの値とサンプル点数Ｋとの比較（ステップＳ５０）が順に行われる。変数ｋの値がサンプル点数Ｋ以下であると判断された場合（ステップＳ５０の判断結果が「ｎｏ」の場合）には、変数ｋの値をインクリメントした後に計算周波数ｆを更新する処理が演算処理装置１４によって再び行われる（ステップＳ２０）。
【００４９】
このように、不安定因子ρが閾値ρｔｈよりも小であり（ステップＳ４３の判断結果が「ｙｅｓ」であり）、且つ変数ｋの値がサンプル点数Ｋ以下である（ステップＳ５０の判断結果が「ｎｏ」である）場合には、順序付け部１４ｂによる順序付け（ステップＳ４１）が行われることなく、計算周波数ｆの更新（ステップＳ２０）、計算周波数ｆに応じた回路方程式の作成（ステップＳ３０）、ＬＵ分解（ステップＳ４２）、及び求解演算（ステップＳ４４）が繰り返し行われる。
【００５０】
以上の処理が行われている最中に、ステップＳ４３で求められる不安定因子ρが閾値ρｔｈ以上になったとすると、ステップＳ４３の判断結果が「ｎｏ」になる。すると、完全ピボッティングによる順序付けが順序付け部１４ｂで行われ（ステップＳ４１）、その順序付けの結果が順序付け結果ＯＲとして補助記憶装置１３に記憶される。その後、補助記憶装置１３に記憶された順序付け結果ＯＲを用いたＬＵ分解がＬＵ分解部１４ｃで行われる（ステップＳ４２）。このように、ステップＳ４３で求められる不安定因子ρが閾値ρｔｈ以上になった場合には、順序付け部１４ｂによる順序付けが行われ、新たに得られた順序付け結果ＯＲを用いてＬＵ分解が行われる。
【００５１】
変数ｋの値がサンプル点数Ｋ以下である場合（ステップＳ５０の判断結果が「ｎｏ」である場合）には、以上説明した一連の処理が繰り返され、サンプル点毎の電力系統の周波数特性が順次求められる。尚、図２に示す一連の処理は、変数ｋの値がサンプル点数Ｋ以下であると判断された場合（ステップＳ５０の判断結果が「ｎｏ」の場合）に終了する。
【００５２】
このように、本実施形態では、計算が終了した周波数と次に計算する周波数との差が小さい場合には、それぞれの周波数で得られる回路方程式の変化も小さいことに着目し、順序付け部１４ｂによる順序付けをサンプル点毎に行わずに、過去に実施した順序付けの結果を再利用するようにしている。このため、周波数特性計算に要する時間を大幅に短縮することが可能である。
【００５３】
また、本実施形態では、順序付け部１４ｂによる順序付けを常に省略するのではなく、必要に応じて順序付け部１４ｂによる順序付けを行うようにしている。具体的には、計算精度の悪化が生じたと評価された場合（不安定因子ρが閾値ρｔｈ以上になった場合）に、順序付け部１４ｂによる順序付けを行うようにしている。このため、計算精度を維持することもできる。
【００５４】
〈計算速度及び計算精度の検証〉
本出願の発明者は、上述した周波数特性解析装置１を用いて電力系統の周波数特性解析を実際に行い、計算速度及び計算精度の検証を行った。本実施形態の周波数特性計算は、計算精度が問題とならない場合に順序付けを省略することによって計算の効率化を行うものであるが、一般に順序付け等の行列演算の計算量は、行列次数Ｎに影響を受ける。このため、以下では、まず行列次数Ｎと計算速度との関係について検証し、次に計算精度及び計算速度ついて従来例と比較しつつ検証する。尚、行列次数Ｎは、ノード数Ｎｎ及びブランチ数Ｎｂにより定まる（Ｎ＝Ｎｎ×２Ｎｂ）ため、回路規模と考えても差し支えない。よって、以下では、Ｎを回路規模と呼ぶ。
【００５５】
（Ａ）回路規模と計算速度との関係
完全ピボッティングによる順序付けの計算量は、回路規模Ｎに影響を受けると考えられるため、回路規模Ｎが変化すると速度向上倍数（従来手法と比較して計算速度が何倍向上したかを示す指標）も変化すると考えられる。そこで、回路規模Ｎと速度向上率との関係を確認することを目的として、回路規模Ｎを変化させた場合の速度向上倍数の算出を行った。
【００５６】
図３は、回路規模と計算速度との関係を検証するために用いた検証用回路を示す図である。図３に示す検証用回路は、発電機Ｇ２〜Ｇ１０、変圧器ＴＲ２〜ＴＲ１０、及び送電線ＴＬ１〜ＴＬ１６を有しており、発電機Ｇ２〜Ｇ１０がそれぞれ変圧器ＴＲ２〜ＴＲ１０を介して送電網（送電線ＴＬ１〜ＴＬ１６からなる送電網）に接続された回路である。この検証用回路の回路規模を図中に示す「ケース１」〜「ケース９」の９通りに変化させ、各々のケースについて周波数特性計算を行った。
【００５７】
例えば、「ケース１」は、発電機Ｇ２、変圧器ＴＲ２、及び送電線ＴＬ１，ＴＬ２からなる規模の回路であり、「ケース２」は、発電機Ｇ２，Ｇ３、変圧器ＴＲ２，ＴＲ３、及び送電線ＴＬ１，ＴＬ２，ＴＬ３，ＴＬ１１からなる規模の回路である。尚、「ケース９」は、図３に示す発電機Ｇ２〜Ｇ１０、変圧器ＴＲ２〜ＴＲ１０、及び送電線ＴＬ１〜ＴＬ１６の全てからなる規模の回路である。「ケース１」〜「ケース９」の各々におけるノード数及びブランチ数を以下の表１に示す。
【００５８】
【表１】
0006650333.tif 000010
【００５９】
図４は、図３中の発電機Ｇ２〜Ｇ１０の模擬に用いた等価回路を示す図である。図４（ａ）に示す通り、発電機Ｇ２〜Ｇ１０の模擬に用いた等価回路は、発電機の起電力Ｅに、次過渡リアクタンスに相当するインダクタンスＬを直列接続した回路である。発電機の次過渡インピーダンスには、周波数が高くなるにつれて表皮効果により損失が大きくなる特性があり、例えば図４（ｂ）に示す等価回路でも模擬できると考えられる。ここでは、計算速度の評価を主目的としていることから、次過渡インピーダンスは損失成分がないと仮定し、簡易的に図４（ａ）に示す等価回路で模擬することとした。発電機Ｇ２〜Ｇ７，Ｇ９のインダクタンスＬは３２．１［μＨ］とし、発電機Ｇ８のインダクタンスＬは６４．２［μＨ］とし、発電機Ｇ１０のインダクタンスＬは１０．７［μＨ］とした。また、電源電圧Ｅは短絡した。
【００６０】
図５は、図３中の変圧器ＴＲ２〜ＴＲ１０の模擬に用いた等価回路を示す図である。図５（ａ）に示す通り、変圧器ＴＲ２〜ＴＲ１０の模擬に用いた等価回路は、Δ−Ｙ接続の理想変圧器に、漏れインダクタンスＬ１が直列に接続された回路である。変圧器の漏れインピーダンスも発電機と同様に表皮効果によって損失が増加する特性があり、例えば図５（ｂ）に示す等価回路でも模擬できると考えられる。ここでは、発電機の模擬と同様の理由で、漏れインピーダンスは損失成分がないと仮定し、簡易的に図５（ａ）に示す等価回路で模擬することとした。
【００６１】
変圧器ＴＲ２〜ＴＲ７，ＴＲ９の漏れインダクタンスＬ１は９．２８［ｍＨ］とし、変圧器ＴＲ８の漏れインダクタンスＬ１は１８．６［ｍＨ］とし、変圧器ＴＲ１０の漏れインダクタンスＬ１は３．１０［ｍＨ］とした。変圧比は、５００［ｋＶ］：２２［ｋＶ］とし、変圧器巻線の対地浮遊静電容量は、何れの変圧器も２０００［ｐＦ］のキャパシタンスで模擬した。送電線ＴＬ１〜ＴＬ１６は、本来であれば周波数依存分布定数線路モデルを用いた現実的な模擬をすべきであるが、ここでの検討では計算速度の評価を主目的としていることから、簡易的に相間の結合を有さないπ型等価回路モデル（１段）を用いた。
【００６２】
図３に示す「ケース１」〜「ケース９」の何れのケースにおいても、図中の母線Ａから系統を見込んだインピーダンスの周波数特性計算を行った。ここで、計算周波数が広範囲であるほど、周波数サンプル間の回路方程式の差が大きくなり、本実施形態における周波数特性計算にとって不利な条件となる。そのため、ここでの検証では、通常実施される周波数範囲よりも広い０．１［Ｈｚ］から１［ＭＨｚ］までとすることで、本実施形態における周波数特性計算の効果を確認することとした。周波数サンプル点は、対数軸上で等間隔に２０００点とした。計算には、現在、一般的に普及しているデスクトップパソコンを用いた。
【００６３】
ここで、一般的な計算機を用いた数値計算では、倍精度演算を用いるのが標準的であり、その有効桁数はおよそ１６桁となる。つまり、前述した不安定因子ρが１６桁を超えた場合には、ＬＵ分解中に実行される計算（前述した（６），（７）式に示す計算）中の差分演算で情報落ちが発生し、計算精度が低下する可能性が生じる。そのため、ここでの計算では、５桁の余裕を見込んで閾値ρｔｈを１０１１とした。尚、ここで用いる「ケース１」〜「ケース９」は、検証用に作成したケースであり、計算結果として得られる周波数特性自体に特別な意味はない点に注意されたい。
【００６４】
図６は、図３に示す「ケース１」〜「ケース９」について周波数特性計算を行った場合に要した計算時間を示す図である。尚、図６中の凡例「従来手法」は、図１２に示す従来の方法で周波数特性計算を行ったときのものであり、凡例「提案手法」は、図２に示す本実施形態の方法で周波数特性計算を行ったときのものである。図６を参照すると、従来手法では、「ケース１」から「ケース９」の順で（つまり、回路規模が大きくなるにつれて）計算時間が増大するのに対し、提案手法では、回路規模が大きくなっても計算時間が殆ど変化しないことが分かる。
【００６５】
図７は、回路規模と速度向上倍数との関係を示す図である。この図７は、図６に示す各々のケースについての計算時間の比（従来手法の計算時間を提案手法の計算時間で除算して得られた値）を速度向上倍数としてプロットしたものである。尚、図６に示す各々のケースの回路規模は、前述した表１から得られる。図７を参照すると、回路規模が大きくなるにつれて速度向上倍数が大きくなることが分かる。特に、図６に示す「ケース９」の場合（回路規模が１３５３の場合）には、速度向上倍数が６０倍近くになっており、速度向上の効果が極めて大きいことが分かる。従って、計算時間が問題となる大規模回路であるほど高速化の効果が大きく、実用的であると考えられる。
【００６６】
図８は、図３に示す「ケース９」について周波数特性計算を行った場合の各サンプル点における不安定因子を示す図である。図８を参照すると、従来手法及び提案手法の不安定因子は、計算周波数が約５０００［Ｈｚ］を超えるあたりまではほぼ同じ値であるが、計算周波数がそれ以上高くなると提案手法では不安定因子が急激に増大することが分かる。また、図８を参照すると、提案手法の不安定因子は、計算周波数が約５００［ｋＨｚ］あたりで急激に減少して、従来手法の不安定因子と同じ値に戻っていることが確認できる。これは、計算周波数が約５００［ｋＨｚ］あたりで提案手法の不安定因子が閾値ρｔｈ（１０１１）を超え、順序付けがやり直されたためである。図８から明らかな通り、全計算過程における初回の計算を含めた順序付けの計算回数は２回であり、全２０００点（周波数サンプル点の数）の計算のうち１９９８回の順序付けが省略されている。
【００６７】
（Ｂ）計算精度及び計算速度の検証
図９は、計算精度及び計算速度を検証するために用いた検証用回路を示す図である。図９に示す検証回路は、発電機Ｇ１〜Ｇ３、変圧器ＴＲ１〜ＴＲ３、周波数変換所ＦＣ、及び送電線ＴＬ１〜ＴＬ３を有しており、発電機Ｇ１〜Ｇ３がそれぞれ変圧器ＴＲ１〜ＴＲ３を介して送電網（送電線ＴＬ１〜ＴＬ３からなる送電網）に接続されるとともに、周波数変換所ＦＣが上記の送電網に接続された２７５［ｋＶ］系統の回路である。この検証用回路は、交直変換器に伴う高調波不安定現象や高調波拡大現象の解析等で行われる交流系統の周波数特性の計算を想定した回路である。尚、図９においては、図示を簡単にするために、単線結線図で描かれているが、実際には全て三相モデルで構成されている。
【００６８】
周波数変換所ＦＣについては、他励式周波数変換所の回路構成（変換器、交流フィルタ、及び調相設備を備える回路構成）を用いた。この周波数変換所ＦＣのパラメータを以下の表２に示す。
【００６９】
【表２】
0006650333.tif 000011
【００７０】
発電機Ｇ１〜Ｇ３は、図４（ｂ）に示す通り、発電機の起電力を表す電圧源に初期過渡インピーダンスが直列に接続された等価回路で模擬した。初期過渡インピーダンスは表皮効果を考慮して図４（ｂ）に示したインダクタンスＬ１，Ｌ２及び抵抗Ｒ２からなる回路で表現した。具体的に、発電機Ｇ１については、インダクタンスＬ１を０．２８３［ｍＨ］とし、インダクタンスＬ２を０．５６６［ｍＨ］とし、抵抗Ｒ２を０．５７７［Ω］とした。また、発電機Ｇ２については、インダクタンスＬ１を０．５６６［ｍＨ］とし、インダクタンスＬ２を１．１３［ｍＨ］とし、抵抗Ｒ２を１．１５［Ω］とした。また、発電機Ｇ３については、インダクタンスＬ１を５４．５［μＨ］とし、インダクタンスＬ２を０．１０９［ｍＨ］とし、抵抗Ｒ２を０．１１１［Ω］とした。高調波の解析であるため、電源電圧Ｅは短絡した。
【００７１】
変圧器ＴＲ１〜ＴＲ３は、図５（ｂ）に示す通り、Δ−Ｙ接続の理想変圧器に、漏れインピーダンスが直列に接続された等価回路で模擬した。変圧器ＴＲ１〜ＴＲ３の漏れインピーダンスは、図５（ｂ）に示したインダクタンスＬ１，Ｌ２及び抵抗Ｒ２からなる回路で表現した。具体的に、変圧器ＴＲ１については、インダクタンスＬ１を２４．２［ｍＨ］とし、インダクタンスＬ２を４．７８［ｍＨ］とし、抵抗Ｒ２を１２．７［Ω］とした。また、変圧器ＴＲ２については、インダクタンスＬ１を４６．５［ｍＨ］とし、インダクタンスＬ２を１１．５［ｍＨ］とし、抵抗Ｒ２を２９．６［Ω］とした。また、変圧器ＴＲ３については、インダクタンスＬ１を４．８４［ｍＨ］とし、インダクタンスＬ２を０．９５５［ｍＨ］とし、抵抗Ｒ２を２．５５［Ω］とした。変圧比は、２７５［ｋＶ］：２２［ｋＶ］とし、変圧器巻線の対地浮遊静電容量は、何れの変圧器も２０００［ｐＦ］のキャパシタンスで模擬した。
【００７２】
送電線ＴＬ１〜ＴＬ３は、全て架空線で、装柱は２７５［ｋＶ］標準装柱であると仮定し、大地抵抗率を１００［Ωｍ］として周波数依存線路モデルで模擬した。また、送電線ＴＬ１〜ＴＬ３の亘長はそれぞれ５［ｋｍ］，１５［ｋｍ］，１００［ｋｍ］で、周波数依存線路モデル作成時の線路定数計算周波数のサンプルは、０．１［Ｈｚ］〜１０［ＭＨｚ］の範囲で４００点とした。
【００７３】
以上の検証用回路について、図９中の点Ｐから見込んだ交流系統のインピーダンスの周波数特性を１［Ｈｚ］から２０００［Ｈｚ］まで１［Ｈｚ］刻みで計算した。計算には、上述した検証（図３に示す検証用回路に対する検証）で用いたデスクトップパソコンを用い、閾値ρｔｈを１０１１とした。
【００７４】
図１０は、図９の検証用回路に対して行った周波数特性計算の計算結果を示す図である。図１０を参照すると、各周波数フィルタの共振次数の周波数においてインピーダンスが低下していることが確認できる。また、図１０では、従来手法の計算結果と提案手法の計算結果との双方を描画しているが、両計算結果は良く一致しており、図１０からはその違いを判別することができない。計算時間は、従来手法が約１５６秒であり、提案手法が約８秒であった。これは、従来手法に比べて約２０倍の計算速度であり、提案手法は従来手法に比べて高速に計算できることが確認できた。
【００７５】
図１１は、図９の検証用回路について周波数特性計算を行った場合の各サンプル点における不安定因子を示す図である。図１１を参照すると、各サンプル点における不安定因子の増加は殆どなく、その推移も波形レベルで一致した。また、不安定因子の大きさは最大で１０５程度であり、提案手法では順序付けは繰り返されず、最も計算量が少ないケースとなった。このように、本検証では、順序付けの結果の再利用が可能であり、それにより計算速度を高速化できることが明らかになった。
【００７６】
以上、本発明の一実施形態による周波数特性解析装置、方法、及びプログラムについて説明したが、本発明は上述した実施形態に制限されることなく、本発明の範囲内で自由に変更が可能である。例えば、上記実施形態では、周波数特性解析装置１が、デスクトップ型のパーソナルコンピュータによって実現される例について説明したが、周波数特性解析装置は、ネットワークに接続されたサーバの形態で実現されても良い。このような形態の場合には、例えばユーザによって操作される端末装置から解析対象及び解析条件がサーバに送信され、サーバで実施された解析の解析結果が端末装置に返信されることとなる。
【００７７】
また、上述した実施形態では、周波数特性解析装置１が周波数特性計算を行う専用の装置として説明したが、周波数特性解析装置１には、周波数特性計算機能以外の他の機能（例えば、瞬時値解析機能）が設けられていても良い。また、上記実施形態では、電力系統の周波数特性計算を行う場合を例に挙げて説明したが、本発明は疎行列演算が必要となる他の分野の解析（例えば、構造解析、流体解析等）にも適用することが可能である。本発明は、計算周波数が変更されると内容が変更されてしまう疎行列に対する疎行列演算を高速化する用途一般に適用することが可能である。
【符号の説明】
【００７８】
１        周波数特性解析装置
１４ａ    回路方程式作成部
１４ｂ    順序付け部
１４ｃ    ＬＵ分解部
１４ｄ    計算精度評価部
１４ｅ    求解演算部
ρ        不安定因子
ρｔｈ      閾値
【図１】
0006650333.tif 000012
【図２】
0006650333.tif 000013
【図３】
0006650333.tif 000014
【図４】
0006650333.tif 000015
【図５】
0006650333.tif 000016
【図６】
0006650333.tif 000017
【図７】
0006650333.tif 000018
【図８】
0006650333.tif 000019
【図９】
0006650333.tif 000020
【図１０】
0006650333.tif 000021
【図１１】
0006650333.tif 000022
【図１２】
0006650333.tif 000023