(57)【要約】【課題】二段対角型陰的ルンゲクッタ法を用いて数値解析を行う際の誤差を高精度に推定し、また、この推定した誤差を用いて計算時間を短縮する。【解決手段】スイッチング素子を含む電子回路の回路素子の特性情報に基づき当該電子回路の回路方程式であるｄｙ／ｄｔ＝ｆ（ｙ）を設定する回路方程式設定処理を実行し、二段対角型陰的ルンゲクッタ法による（ａ）式及び（ｂ）式を用いて前記回路方程式を解き、第ｎステップ目（ｎは自然数）の解ｙｎを算出する解析処理を実行し、（ｄ）式の係数ａ１，ａ２を求め、（ｅ）式の係数ｂ１，ｂ２，ｂ３を求め、解ｙｎ−１のテイラー展開を用いて求めた時刻ｔｎにおける解析解ｙ＊ｎと（ｅ）式のｙｎとの差に前記係数ｂ１，ｂ２，ｂ３を代入することで解ｙｎの誤差εを推定する（ｆ）式を算出し、当該（ｆ）式により誤差εを推定する誤差推定処理を実行する。【選択図】図３
2010282359.tif 000002
【特許請求の範囲】
【請求項１】
回路方程式設定手段と、解析手段と、誤差推定手段とを備えたコンピュータにおいて、  前記回路方程式設定手段が、スイッチング素子を含む電子回路の回路素子の特性情報に基づき当該電子回路の回路方程式であるｄｙ／ｄｔ＝ｆ（ｙ）を設定する回路方程式設定処理を実行し、
前記解析手段が、二段対角型陰的ルンゲクッタ法による以下の（ａ）式及び（ｂ）式を用いて前記回路方程式を解き、第ｎステップ目（ｎは自然数）の解ｙｎを算出する解析処理を実行し、
前記誤差推定手段が、（ｂ）式の右辺に（ｄ）式を代入し、テイラー展開することにより得られる式と（ｄ）式とを比較することにより、（ｄ）式の係数ａ１，ａ２を求め、この（ａ）式及び（ｂ）式より導かれる（ｃ）式の右辺に（ｄ）式、（ｅ）式及び前記係数ａ１，ａ２を代入し、テイラー展開することにより得られる式と（ｅ）式を比較することにより、（ｅ）式の係数ｂ１，ｂ２，ｂ３を求め、解ｙｎ−１のテイラー展開を用いて求めた時刻ｔｎにおける解析解ｙ＊ｎと（ｅ）式のｙｎとの差に前記係数ｂ１，ｂ２，ｂ３を代入することで解ｙｎの誤差εを推定する（ｆ）式を算出し、当該（ｆ）式により誤差εを推定する誤差推定処理を実行する
ことを特徴とする回路解析方法。
【数１】
2010282359.tif 000036
【数２】
2010282359.tif 000037
【数３】
2010282359.tif 000038
【数４】
2010282359.tif 000039
【数５】
2010282359.tif 000040
【数６】
2010282359.tif 000041
【請求項２】
請求項１に記載する回路解析方法において、
前記（ｆ）式のｄ２ｆ／ｄｔ２は
【数７】
2010282359.tif 000042
であることを特徴とする回路解析方法。
【請求項３】
請求項１又は請求項２に記載する回路解析方法において、
前記解析処理では、前記スイッチング素子の状態変化と、前記誤差推定処理により算出した解ｙｎ−１の誤差εに基づいて時間刻み幅ｈを決定する時間刻み幅設定処理を実行した後、解ｙｎを算出する
ことを特徴とする回路解析方法。
【請求項４】
請求項１〜請求項３の何れか一項に記載する回路解析方法をコンピュータに実行させるための回路解析プログラム。
【請求項５】
請求項４に記載する回路解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項６】
請求項１〜請求項３の何れか一項に記載する回路解析方法を実行するように構成されたコンピュータを有することを特徴とする回路解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回路解析方法、回路解析装置、回路解析プログラム、回路解析プログラムを記録した記録媒体に関する。
【背景技術】
【０００２】
近年、電気エネルギーの利用において、電圧、電流、周波数などの変換を行う半導体電力変換装置は、民生、産業、運輸分野の様々な機器に応用され、その高性能化を実現する上で重要な技術となっている。一方、半導体電力変換装置の開発現場においては、多大なコストを要する模型試験の回数を削減したり、半導体電力変換装置の挙動を解析したりするなどの手段として、数値解析（シミュレーション）は不可欠なものとなっている。
【０００３】
特に、トランジスタ等のスイッチング素子を含む電子回路の過渡現象解析に際しては、スイッチング素子の特性を詳細に考慮し、シミュレーションを行う必要がある。一般に高精度なシミュレーションを行うと、計算負荷が大きくなり計算時間が増大するという問題がある。
【０００４】
そこで、局所打ち切り誤差を推定し、この誤差が所定値以下となるように、シミュレーションを実行する際の時間刻み幅を可変とする方法がある（例えば、特許文献１参照）。例えば、スイッチング動作に掛かる時間は数十ナノ秒程度であるので、その時間においては時間刻み幅を数ナノ秒程度にし、それ以外の時間においては１００ナノ秒〜１マイクロ秒程度とすることで、全体として計算時間を短縮することができる。
【０００５】
ここで、シミュレーションに用いられる数値積分手法である二段対角型陰的ルンゲクッタ法（以下、2S-DIRK（2-Stage Diagonally Implicit Runge-Kutta）とも表記する。）は、台形則と同じ二次の精度を持ちながら、台形則のような数値振動の持続が生じないことが知られている（例えば、非特許文献１、非特許文献２参照。）。したがって、2S-DIRKは、スイッチング素子を含む電子回路の過渡現象をシミュレーションする際に用る数値積分手法として適していると考えられるが、その打ち切り誤差を高精度に推定する方法は確立されていないのが現状である。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００７−３２８３８６号公報（請求項６）
【非特許文献】
【０００７】
【非特許文献１】Taku Noda, Kiyoshi Takenaka, Toshio Inoue, “Numerical Integration by the 2-Stage Diagonally Implicit Runge-Kutta Method for Electromagnetic Transient Simulations”, IEEE Transactions Power Delivery, vol. 24, No. 1, pp. 390-399, 2009
【非特許文献２】井上俊夫、田中和幸、熊野照久「電力系統長時間動特性シミュレーション手法−系統電圧特性に着目した基本プログラムの開発−」、電力中央研究所  研究報告T01051, 2002
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明は、上記従来技術に鑑み、二段対角型陰的ルンゲクッタ法を用いて数値解析を行う際の誤差を高精度に推定し、また、この推定した誤差を用いて計算時間を短縮することができる回路解析方法、回路解析装置、回路解析プログラム、回路解析プログラムを記録した記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記目的を達成するための本発明の第１の態様は、回路方程式設定手段と、解析手段と、誤差推定手段とを備えたコンピュータにおいて、前記回路方程式設定手段が、スイッチング素子を含む電子回路の回路素子の特性情報に基づき当該電子回路の回路方程式であるｄｙ／ｄｔ＝ｆ（ｙ）を設定する回路方程式設定処理を実行し、前記解析手段が、二段対角型陰的ルンゲクッタ法による以下の（ａ）式及び（ｂ）式を用いて前記回路方程式を解き、第ｎステップ目（ｎは自然数）の解ｙｎを算出する解析処理を実行し、前記誤差推定手段が、（ｂ）式の右辺に（ｄ）式を代入し、テイラー展開することにより得られる式と（ｄ）式とを比較することにより、（ｄ）式の係数ａ１，ａ２を求め、この（ａ）式及び（ｂ）式より導かれる（ｃ）式の右辺に（ｄ）式、（ｅ）式及び前記係数ａ１，ａ２を代入し、テイラー展開することにより得られる式と（ｅ）式を比較することにより、（ｅ）式の係数ｂ１，ｂ２，ｂ３を求め、解ｙｎ−１のテイラー展開を用いて求めた時刻ｔｎにおける解析解ｙ＊ｎと（ｅ）式のｙｎとの差に前記係数ｂ１，ｂ２，ｂ３を代入することで解ｙｎの誤差εを推定する（ｆ）式を算出し、当該（ｆ）式により誤差εを推定する誤差推定処理を実行することを特徴とする回路解析方法にある。
【００１０】
【数１】
2010282359.tif 000003
【数２】
2010282359.tif 000004
【数３】
2010282359.tif 000005
【数４】
2010282359.tif 000006
【数５】
2010282359.tif 000007
【数６】
2010282359.tif 000008
【００１１】
かかる第１の態様では、二段対角型陰的ルンゲクッタ法を用いて誤差の小さい高精度な数値解析を実現することができる。また、誤差を小さくできるため、その誤差を基準にして時間刻み幅ｈを決定する数値解析を行う際に、計算時間を短縮することができる。
【００１２】
本発明の第２の態様は、第１の態様に記載する回路解析方法において、前記（ｆ）式のｄ２ｆ／ｄｔ２は
【数７】
2010282359.tif 000009
であることを特徴とする回路解析方法にある。
【００１３】
かかる第２の態様では、４つの解に基づいて誤差を推定することができる。
【００１４】
本発明の第３の態様は、第１又は第２の態様に記載する回路解析方法において、前記解析処理では、前記スイッチング素子の状態変化と、前記誤差推定処理により算出した解ｙｎ−１の誤差εに基づいて時間刻み幅ｈを決定する時間刻み幅設定処理を実行した後、解ｙｎを算出することを特徴とする回路解析方法にある。
【００１５】
かかる第３の態様では、電子回路を過渡状態にし得るスイッチング素子の状態変化に着目して時間刻み幅ｈを決定すること、電子回路中の他の素子を考慮に入れた場合よりも計算時間を短縮することができる。
【００１６】
本発明の第４の態様は、第１〜第３の何れか一つの態様に記載する回路解析方法をコンピュータに実行させるための回路解析プログラムにある。
【００１７】
本発明の第５の態様は、請求項４に記載する回路解析プログラムを記録したコンピュータ読み取り可能な記録媒体にある。
【００１８】
本発明の第６の態様は、請求項１〜請求項３の何れか一項に記載する回路解析方法を実行するように構成されたコンピュータを有することを特徴とする回路解析装置にある。
【００１９】
かかる第４〜第６の態様は、いずれも、これらを用いて上記回路解析を実施することができ、二段対角型陰的ルンゲクッタ法を用いて高精度な誤差の推定を実現するとともに、計算時間を短縮することができる。
【発明の効果】
【００２０】
本発明によれば、二段対角型陰的ルンゲクッタ法を用いて誤差の小さい高精度な数値解析を実現するとともに、当該誤差に基づいた可変時間刻み幅を採用することで計算時間を短縮することができる回路解析方法、回路解析装置、回路解析プログラム、回路解析プログラムを記録した記録媒体が提供される。
【図面の簡単な説明】
【００２１】
【図１】本発明に係る回路解析装置の構成を示すブロック図である。
【図２】本発明に係る回路解析装置のコンピュータのハードウェア構成を示すブロック図である。
【図３】本発明に係る回路解析方法のフローを示す図である。
【図４】本発明に係る回路解析方法の誤差推定処理に関するフローを示す図である。
【図５】比較例に係る誤差推定方法のイメージである。
【図６】本発明に係る回路解析方法の解析対象であるＲＣ回路図である。
【図７】本発明に係る回路解析方法の誤差と比較例（推定手法２）の誤差とを示すグラフである。
【図８】本発明に係る回路解析方法の解析対象であるＲＬＣ回路図である。
【図９】本発明に係る回路解析方法の誤差と比較例（推定手法２）の誤差とを示すグラフである。
【図１０】本発明に係る回路解析方法の解析対象であるチョッパ回路図である。
【図１１】本発明に係る回路解析方法によるＭＯＳＦＥＴの電圧を示すグラフである。
【図１２】本発明に係る回路解析方法による時間刻み幅の推移を示すグラフである。
【図１３】本発明に係る回路解析方法の解析対象であるインバータ回路図である。
【図１４】本発明に係る回路解析方法によるインバータ出力電流及びＩＧＢＴ電圧を示すグラフである。
【発明を実施するための形態】
【００２２】
以下、本発明の実施形態に係る回路解析方法及び回路解析装置について説明する。本実施形態に係る回路解析装置１は、図１に示すように、回路方程式設定処理、解析処理及び誤差推定処理を実行する各手段１１〜１３を備えたコンピュータ２０を有する装置である。図２に示すように、コンピュータ２０は、ＣＰＵ２１、ＲＡＭ２２、ＲＯＭ２３、ハードディスク２４等の記憶手段を備え、キーボードや記憶媒体読み取り装置などの入力装置３０と、ディスプレイ、プリンタ等の出力装置４０が接続されている。
【００２３】
図１に示すように、コンピュータ２０が実行可能な手段は、後述する回路解析方法における各処理を実行する手段であり、具体的には、回路方程式設定手段１１と、解析手段１２と、誤差推定手段１３とを有している。なお、これらの各手段１１〜１３は、後述する回路解析方法をコンピュータ２０に実行させるための回路解析プログラム１０として作成したものを読み込んで実行可能としたものである。
【００２４】
ここで、図３を用いて、回路解析方法の処理を詳細に説明する。図示するように、まず、スイッチング素子を含む電子回路の回路素子の特性情報に基づき当該電子回路の回路方程式であるｄｙ／ｄｔ＝ｆ（ｙ）を設定する回路方程式設定処理を実行する（ステップＳ１）。この処理は、回路方程式設定手段１１により実行される。具体的には、ハードディスク２４に記録された回路データ、又は、入力装置３０を介して解析対象の回路の回路データをＲＡＭ２２に読み込み、この回路データに基づいて回路方程式を作成する。
【００２５】
回路データは、電子回路を構成する各電子部品の特性とその接続関係を表す情報等を含む。回路方程式は、解析対象回路の接点電位、素子電流及び素子電圧を変数とする連立方程式であり、解析対象回路の接点の数をＮｓ、素子の数をＮｂとすると、Ｎｓ＋２Ｎｂ個の数の方程式を有する。
【００２６】
次に、二段対角型陰的ルンゲクッタ法により、回路方程式を解く解析処理を行う。具体的には、第ｎステップ目（ｎは自然数）の解を解ｙｎとすると、解析手段１２は、初期解ｙ０等の回路方程式を解くために必要な初期値を設定する初期値設定処理（ステップＳ２）、第ｎステップ目の解ｙｎの計算処理で採用する時間刻み幅ｈの決定処理（ステップＳ３）、時刻ｔｎにおける解ｙｎの計算処理（ステップＳ４）を順次実行する。
【００２７】
そして、誤差推定手段１３が、解ｙｎについての誤差を推定する誤差推定処理を実行する（ステップＳ５）。最後に、計算終了時刻を越えたか否かを判別し（ステップＳ６）、すなわち、時刻ｔｎが所定の時刻を越えていないならば（ステップＳ６：Ｎｏ）、第ｎ＋１ステップ目として、ステップＳ３〜Ｓ５を実行する。一方、時刻ｔｎが所定の時刻を越えているならば（ステップＳ６：Ｙｅｓ）、計算を終了する。
【００２８】
以下、順番が前後するが、まず時刻ｔｎにおける解ｙｎの計算処理について詳細に説明する。すなわち、時間刻み幅ｈの決定処理（ステップＳ３）により時間刻み幅ｈが決定されたものとして、当該計算処理について説明する。
【００２９】
本発明に係る回路解析方法では、2S-DIRKにより回路方程式を解く。このとき、時間刻み幅ｈで、時刻ｔｎ−１からｔｎに進むための数式は次の通りとなる。
【００３０】
まず、数９の式を用いて、一時的な解
【数８】
2010282359.tif 000010
を算出する。なお、「ｓ」は、2S-DIRKでは、ｓ＝１−１／SQRT(2)（SQRTは平方根を求める関数）であり、前述の初期値設定処理（ステップＳ２）で予め入力し、ＲＡＭ２２に記憶しておく。
【００３１】
【数９】
2010282359.tif 000011
【００３２】
次に、数１０の式を用いて、時刻ｔｎにおける解ｙｎを求める。2S-DIRKでは、α＝−SQRT(2)、β＝１＋SQRT(2)であり、前述の初期値設定処理（ステップＳ２）で予め入力し、ＲＡＭ２２に記憶しておく。
【００３３】
【数１０】
2010282359.tif 000012
【００３４】
なお、数９及び数１０に記載の各式から解ｙｎは数１１の（ｃ）式のように表すこともできる。
【００３５】
【数１１】
2010282359.tif 000013
【００３６】
次に、上記計算により求めた解ｙｎについての誤差εを推定する誤差推定処理（ステップＳ５）について説明する。
【００３７】
まず、解
【数１２】
2010282359.tif 000014
を（ｄ）式のように表し、（ｂ）式の右辺に（ｄ）式を代入し、テイラー展開することにより得られる式と（ｄ）式を比較することにより（ｄ）式の係数ａ１，ａ２を求める。
【００３８】
【数１３】
2010282359.tif 000015
【００３９】
ここで、（ｄ）式は次の（ｄ１）式のように表せる。
【数１４】
2010282359.tif 000016
【００４０】
（ｂ）式の右辺の解に（ｄ１）式を代入すると、以下の（ｄ２）式が得られる。
【数１５】
2010282359.tif 000017
【００４１】
（ｄ２）式を解ｙｎ−１でテイラー展開すると、次の（ｄ３）式が得られる。
【数１６】
2010282359.tif 000018
【００４２】
（ｄ１）式の右辺と（ｄ３）式の右辺を比較すると以下の（ｄ４）式が得られ、ａ１＝１，ａ２＝１と決定できる。
【数１７】
2010282359.tif 000019
【００４３】
次に、解ｙｎを（ｅ）式のように表す。
【数１８】
2010282359.tif 000020
【００４４】
（ｅ）式は、数１４と同様にして（ｅ１）式のように表すことができる。
【数１９】
2010282359.tif 000021
【００４５】
ここで、（ｃ）式に（ｄ１）式及び（ｅ１）式を代入すると、（ｅ２）式が得られる。
【数２０】
2010282359.tif 000022
【００４６】
（ｅ２）式を解ｙｎ−１でテイラー展開すると、（ｅ３）式が得られる。
【数２１】
2010282359.tif 000023
【００４７】
ここで、物理現象に対して、時間刻み幅ｈを十分細かく取っていて線形と見なせる場合、
【数２２】
2010282359.tif 000024
とすることができるので、（ｅ１）式の右辺と、（ｅ３）式の右辺とを比較すると（ｅ４）式が得られる。
【００４８】
【数２３】
2010282359.tif 000025
【００４９】
よって、係数ｂ１＝１、ｂ２＝１／２、ｂ３≒０．２０７１と決定できる。ここで、時刻ｔｎでの解析解ｙ＊ｎは、ｙｎ−１のテイラー展開を用いて（ｈ）式で表せる。
【数２４】
2010282359.tif 000026
【００５０】
（ｅ３）式にｂ１＝１、ｂ２＝１／２、ｂ３≒０．２０７１を代入した式と、（ｈ）式との差により誤差εを推定する（ｆ）式が求まる。
【数２５】
2010282359.tif 000027
【００５１】
なお、（ｆ）式中の、ｄ２ｆ／ｄｔ２については、次のようにして求めることができる。ｙｎ−３、ｙｎ−２、ｙｎ−１をｙｎに関してテイラー展開すると、数２６の通りとなる。
【数２６】
2010282359.tif 000028
【００５２】
数２６の４式から、
【数２７】
2010282359.tif 000029
となり、
【数２８】
2010282359.tif 000030
と表すことができる。
【００５３】
以上のようにして、時刻ｔｎにおける解ｙｎとその誤差εを計算した。その後は、前述したように、計算終了時刻を越えたか否かにより、計算を終了するか、時間刻み幅ｈの決定処理を行う（図３参照）。
【００５４】
時間刻み幅ｈは、スイッチング素子の状態変化、例えば、ゲート電圧の変化と、誤差εに基づいて定める。図４を用いて、この時間刻み幅ｈの決定処理について詳細に説明する。まず、初期値設定処理において、時間刻み幅の最小値である最小刻み幅ｈｍｉｎ、最大値である最大刻み幅ｈｍａｘ、刻み延長禁止ステップ数（Ｓａ、Ｓｂ）、許容できる誤差の最大値である許容誤差ｒを予め入力しＲＡＭ２２に記憶しておく。
【００５５】
次に、ゲート電圧源や、回路中のスイッチング素子に状態変化があるか否かを確認する（ステップＳ１０）。
【００５６】
ここで、スイッチング素子の状態変化に着目するのは、単に、誤差εのみに基づいて時間刻み幅ｈを決定しようとすると、計算時間が膨大となるからである。詳言すると、特に、回路素子は数百から数千個存在するような複雑な回路においては、回路素子にはインダクタやキャパシタといった動的素子の数も相当数になるが、これら全てについて解を保存し（（ｇ）式にも示したように４点の解を必要とする）、各素子についての誤差全てを対象として時間刻み幅ｈを決定することになるからである。
【００５７】
一方、スイッチング素子を含む電子回路、特にパワーエレクトロニクス回路の特徴は、スイッチング素子のON／OFFによる転流動作を行うことにより電力変換を行っていることや、パワーエレクトロニクス回路において、一般的に最も時定数の短い過渡現象は、スイッチング素子のスイッチングによる電流の急変であることが挙げられる。
【００５８】
よって、パワーエレクトロニクス回路をはじめとするスイッチング素子を含む電子回路においては、スイッチング素子の状態、例えば電流の変化を監視していれば回路が過渡状体にあるのか、過渡状体収束後の安定状態にあるかが判定できる。
【００５９】
このような理由により、時間刻み幅ｈを決定する条件の一つとしてスイッチング素子の状態変化を用いることで、計算時間の短縮を図ることができる。
【００６０】
スイッチング素子の状態、例えばゲート電圧の状態変化があった場合（ステップＳ１０：Ｙｅｓ）、時間刻み幅ｈをｈｍｉｎとし、現在のステップからｓａステップ後までの間は、ｈ＝ｈｍｉｎとする（ステップＳ１１）。
【００６１】
次に、状態変化がない場合（ステップＳ１０：Ｎｏ）であって、現在のステップが、時間刻み幅変更禁止時刻に無い場合、誤差εｍａｘと許容誤差ｒを比較し、εｍａｘ＜０．１×ｒかつｈ＜ｈｍａｘならば（ステップＳ１２：Ｙｅｓ）、時間刻み幅ｈを２倍にする（ステップＳ１３）。なお、誤差εｍａｘとは、直近の数ステップ、例えば１０ステップにおける誤差εの最大値である。このような直近１０ステップの誤差εの最大値である誤差εｍａｘを用いるのは、解析対象になっている物理現象の時定数と比べて、時間刻み幅ｈが粗いなどの理由により、或るステップで誤差を過小に評価してしまった場合に、時間刻み幅ｈを延長するのを防止するためである。ステップＳ１２の条件が成立しない場合（ステップＳ１２：Ｎｏ）、εｍａｘ＞ｒかつｈ＞ｈｍｉｎならば（ステップＳ１４：Ｙｅｓ）、時間刻み幅ｈを半分にする（ステップＳ１５）。なお、ステップＳ１３又はステップＳ１５で時間刻み幅ｈを変更した場合は、その後、ｓｂステップ経過するまでは時間延長を禁止する。
【００６２】
なお、ステップＳ１２において、εｍａｘ＜０．１×ｒとしたのは、（ｆ）式より、2S-DIRKにおいて、時間刻み幅ｈが２倍になれば誤差が８倍となるからである。つまり、現在の誤差推定値が許容誤差ｒの１／８以下になっていれば、時間刻み幅ｈを２倍にしても次のステップでの誤差が許容誤差ｒを超えることはないと考えられるからである。上記の例では、さらに余裕を見て、誤差εが許容誤差の１／１０以下になっていれば時間刻み幅ｈを２倍にすることとした。
【００６３】
以上に説明した回路解析方法及び当該方法を実行するように構成された回路解析装置によれば、スイッチング素子を含む電子回路の回路方程式を2S-DIRKを用いて解くとともに、その誤差を推定することができる。以下、本発明に係る回路解析方法により推定された誤差の精度と、計算時間が短縮されたことを比較例との比較により説明する。また、比較に際しては、近似的に求めた誤差の真値との比較も行う。
【００６４】
［比較例］
まず、本発明に係る回路解析方法との比較対象となる誤差の推定方法を説明する。2S-DIRKでは、時刻ｔｎ−１からｔｎに進む際、両者の中間的な時刻
【００６５】
【数２９】
2010282359.tif 000031
での解
【数３０】
2010282359.tif 000032
を求め（（ａ）式）、その後、時刻ｔｎでの解ｙｎを求める。また、解ｙｎ−１から解
【数３１】
2010282359.tif 000033
を求める際に使用する計算手法は後退オイラー法と同一の手法である。よって、近似的には、図５に示すように、解ｙｎ−１と解
【数３２】
2010282359.tif 000034
を結ぶ直線より求める時刻ｔｎでの解ｙｎｂは、解ｙｎ−１から解ｙｎを時間刻み幅ｈの後退オイラー法によって求めた解の近似解と考えることもできる。解ｙｎと解ｙｎｂとの差をεｂとする。解析対象となっている現象に対して時間刻み幅ｈが十分に小さければ現象は線形近似でき、解ｙｎと解ｙｎｂの差は小さくなる。後退オイラー法は一次の精度であり、2S-DIRKは二次の精度であるので、ｙｎｂと比較してｙｎは十分真値に近い。よって、εｂは後退オイラー法の局所打ち切り誤差の近似値を求めていることとなる。εｂは直接2S-DIRKの局所打ち切り誤差を求めるものではないものの、2S-DIRKの局所打ち切り誤差は通常εｂ以下となっていると考えられるので、εｂは局所打ち切り誤差の目安にできる。
【００６６】
［誤差の真値εｔ］
近似的に誤差の真値εｔを次のステップで算出する。
（１）時間刻み幅ｈで１回積分を行うことにより、時刻ｔｎ−１から計算値ｙｎを求める。
（２）時間刻み幅０．０１ｈで１００回積分を行うことにより、時刻ｔｎ−１から近似的な真値ｙ＊ｎを求める。ｙ＊ｎは厳密には真値ではないが、時間刻み幅ｈを上記（１）で求める場合の１／１００としているので、（１）で求めた計算値ｙｎと比べて十分に真値に近い。
（３）（１）、（２）で求めたｙｎとｙ＊ｎとの差により、１ステップの積分で生じた打ち切り誤差の真値εｔを求める。
（４）時刻ｔｎでの解をｙｎとして、同様に時刻ｔｎ＋１での解ｙｎ＋１を算出する。
【００６７】
［誤差の推定精度の比較１］
図６は、解析対象の電子回路の一例であるＲＣ回路である。図７は、このＲＣ回路について計算した、本発明に係る回路解析方法により誤差εと、比較例による誤差εｂと、誤差の真値εｔとを比較したものである。図中、「推定手法１」は、本発明に係る回路解析方法による誤差、「推定手法２」は、比較例による誤差εｂを示している（以下、他の比較図においても同様である。）。詳細には、図７（ａ）は時間刻み幅ｈを１ｍｓとし、図７（ｂ）は時間刻み幅ｈを０．１ｓとした場合の誤差を比較したものである。
【００６８】
図７（ａ）においては、推定手法１の誤差εは、完全に誤差の真値εｔと重なっており、図７（ｂ）においては、推定手法１の誤差εは、電圧印加直後（ｔ＝０付近）において誤差が見られるものの、その後はほぼ誤差の真値εｔと一致している。また、いずれの時間刻み幅ｈであっても、推定手法２の誤差εｂは、推定手法１及び真値εｔと比べて打ち切り誤差を過大に見積もっている。また、図７（ａ）と図７（ｂ）とを比較すると分かるように、時間刻み幅ｈを大きくすると、推定手法１による誤差εは、誤差の真値εｔよりも過大に誤差を見積もっている。これは、誤差を推定するに当たり、安全側に見積もっていることとなる。
【００６９】
［誤差の推定精度の比較２］
図８は、解析対象の電子回路の一例であるＲＬＣ回路である。図９は、このＲＬＣ回路について計算した、本発明に係る回路解析方法により誤差εと、比較例による誤差εｂと、誤差の真値εｔとを比較したものである。詳細には、図９（ａ）は時間刻み幅ｈを１０μｓとし、図９（ｂ）は時間刻み幅ｈを２００μｓとした場合の誤差を比較したものである。
【００７０】
前記した「誤差の推定精度の比較１」と同様に、図９（ａ）においては、推定手法１の誤差εは、完全に誤差の真値εｔと重なっており、図９（ｂ）においては、推定手法１の誤差εは、ほぼ誤差の真値εｔと一致している。また、いずれの時間刻み幅ｈであっても、推定手法２の誤差εｂは、推定手法１及び真値εｔと比べて打ち切り誤差を過大に見積もっている。
【００７１】
ここで、本発明に係る回路解析方法（推定手法１）と推定手法２との特徴を表１に示す。
【００７２】
【表１】
2010282359.tif 000035
【００７３】
[計算時間の比較]
前記したように、本発明に係る回路解析方法によれば、推定手法２よりも高精度に誤差を推定することができる。
【００７４】
このことは、より適切に時間刻み幅ｈを変更して計算時間を短縮することに寄与する。なぜならば、（１）スイッチング素子を含む電子回路、特に、パワーエレクトロニクス回路では、非常に短い周期でスイッチングを行っているため、一度に時間刻み幅ｈを延長できる区間が短い、（２）実際には時間刻み幅ｈが延長できるにも関わらず、推定手法２による打ち切り誤差の推定は誤差を過大に見積もるため、延長できない（誤差εｂが許容誤差ｒに近くなるためである。図４ステップ１３参照）、という理由によるからである。
【００７５】
このため、推定手法２によれば、回路が過渡状態収束後の安定状態になっても、時間刻み幅ｈを大きくできないことにより計算時間が短縮できない可能性が高い。一方、本発明に係る回路解析方法によれば、誤差を高精度に推定できるため、過渡状態収束後の安定状態に追随して時間刻み幅ｈを短縮することができるので、計算時間を短縮することができる。
【００７６】
以下、本発明に係る回路解析方法による解を求めるに際し、固定時間刻みで行った場合とを比較し、可変時間刻みで行った場合においても高精度に解を計算できることを示す。
【００７７】
[固定時間刻みとの比較１]
図１０は、本比較における解析対象のチョッパ回路である。同図の丸印部分は、スイッチング素子（ＭＯＳＦＥＴ）の電流を監視する部分を示している。図１１は、スイッチング素子の電圧−時間のグラフである。同図には、可変時間刻み（本発明の回路解析方法）と固定時間刻みのそれぞれにより計算された電圧値が示されている。
【００７８】
シミュレーション条件は、ｈｍｉｎ＝２ｎｓ、ｈｍａｘ＝１２８ｎｓ、ｒ＝０．０１Ａ、ｓａ＝１００、ｓｂ＝１２である。最小時間刻みを２ｎｓに設定したのは、ＭＯＳＦＥＴのスイッチングにかかる期間は、通常数十ｎｓなので、ＭＯＳＦＥＴのスイッチング時の波形をシミュレーションで正確に再現しようという場合、時間刻み幅ｈは数ｎｓ程度に設定する必要があるからである。
【００７９】
図１１（ａ）に示すように、可変時間刻みによる計算値は、固定時間刻みと区別がつかないほど一致している。また、図１１（ａ）を拡大した図１１（ｂ）においては、可変時間刻みとしたことにより、ゲート電圧の立ち上がりが遅れたため、波形は完全には一致していないが、ほぼ同様の波形となっている。
【００８０】
また、可変時間刻みで計算した時間刻み幅の推移を図１２に示す。同図に示すように、時間刻み幅ｈが細かくなっているのは、ゲートドライブ電圧が変化し、回路中のＭＯＳＦＥＴがスイッチングを行っているときであり、スイッチング期間が終了し、回路の過渡現象が収まるにつれて時間刻み幅ｈも延長されていることが示されている。
【００８１】
なお、本比較では、固定時間刻みで１８４秒、可変時間刻み（本発明の回路解析方法）では５秒の時間を要した。すなわち、本回路において本発明に係る回路解析方法では、３７倍高速化されたこととなる。計算にはＣＰＵがＩｎｔｅｌ（Ｒ）Ｃｏｒｅ（ＴＭ）２ＱｕａｄＣＰＵ、２．４０ＧＨｚのコンピュータを用いた。
【００８２】
[固定時間刻みとの比較２]
本比較では、スイッチングを行う半導体素子数が多く、スイッチング周波数が高い回路で計算時間の短縮効果を検証する。図１３は、本比較における解析対象のインバータ回路である。同回路は、昇圧チョッパと単相インバータを組み合わせたもので、昇圧チョッパにてＭＯＳＦＥＴを１つ、単相インバータにてＩＧＢＴを４つと、合計５個の素子を用いてスイッチングを行っている。また、昇圧チョッパのキャリア周波数は７２ｋＨｚ、単相インバータのキャリア周波数は１８ｋＨｚと、キャリア周波数が高く設定されている。なお、同図の丸印部分は、電流の打ち切り誤差を監視する部分である。
【００８３】
図１４（ａ）に示すように、本発明に係る回路解析方法によるインバータ出力電流は、固定時間刻みの場合と比較して、区別がつかないほど一致している。また、図１４（ｂ）、（ｃ）に示すように、ＩＧＢＴ電圧の波形についても、本発明に係る回路解析方法と固定時間刻みの双方とも、区別がつかないほど一致している。
【００８４】
０ｍｓから３０ｍｓまでｈｍｉｎ＝２ｎｓ、ｈｍａｘ＝１２８ｎｓ、ｒ＝０．０１Ａ、ｓａ＝２５０、ｓｂ＝１２と設定して計算を行った場合、固定時間刻みで１８３１３秒、可変時間刻みで５２４９秒の時間を要した。すなわち、本回路において本発明に係る回路解析方法では、３．４８倍高速化されたこととなる。
【００８５】
以上に説明したように、本発明に係る回路解析方法、及び同方法を実行するように構成された回路解析装置では、二段対角型陰的ルンゲクッタ法を用いて誤差の小さい高精度な数値解析を実現するとともに、当該誤差に基づいた可変時間刻み幅を採用することで計算時間を短縮することができる。
【産業上の利用可能性】
【００８６】
本発明は、スイッチング素子を含む電子回路の過渡現象解析を行う産業分野で利用することができる。
【符号の説明】
【００８７】
１            回路解析装置
１０          回路解析プログラム
１１          回路方程式設定手段
１２          解析手段
１３          誤差推定手段
２０          コンピュータ
２１          ＣＰＵ
２２          ＲＡＭ
２３          ＲＯＭ
２４          ハードディスク
３０          入力装置
４０          出力装置
【図１】
2010282359.tif 000043
【図２】
2010282359.tif 000044
【図３】
2010282359.tif 000045
【図４】
2010282359.tif 000046
【図５】
2010282359.tif 000047
【図６】
2010282359.tif 000048
【図８】
2010282359.tif 000049
【図７】
2010282359.tif 000050
【図９】
2010282359.tif 000051
【図１０】
2010282359.tif 000052
【図１１】
2010282359.tif 000053
【図１２】
2010282359.tif 000054
【図１３】
2010282359.tif 000055
【図１４】
2010282359.tif 000056