(57)【要約】【課題】  陰的ルンゲクッタ法を利用してシミュレーション解析を行う際のノイズとなるスパイクの発生を防止し得る非線形要素を含む系の解析装置を提供する。【解決手段】  陰的ルンゲクッタ法により非線形要素を有する系の解析を行う非線形要素を含む系の解析装置であって、動作点が非線形要素の特性曲線上に乗るように、１段目の計算のみ前記非線形要素の特性を所定量に縮小した特性を用いて計算を行う特性補償手段３Ａを有する。【選択図】  図１０
2014222486.tif 000002
【特許請求の範囲】
【請求項１】
陰的ルンゲクッタ法により非線形要素を有する系の解析を行う非線形要素を含む系の解析装置であって、
動作点を前記非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更した特性を用いて計算を行う特性補償手段を有することを特徴とする非線形要素を含む系の解析装置。
【請求項２】
請求項１に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法において、変数変換により、各段における移動距離をＬ倍に拡大又は縮小する操作に対し、前記非線形要素の特性曲線を１／Ｌとすることにより動作点が非線形要素の特性曲線上に乗るように前記非線形要素の特性を変更したことを特徴とする非線形要素を含む系の解析装置。
【請求項３】
請求項１または請求項２に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法は二段対角型陰的ルンゲクッタ法であることを特徴とする非線形要素を含む系の解析装置。
【請求項４】
請求項１〜請求項３の何れか一つに記載する非線形要素を含む系の解析装置において、
前記系は電気回路であり、前記非線形要素は前記電気回路を構成する非線形素子であることを特徴とする非線形要素を含む系の解析装置。
【請求項５】
陰的ルンゲクッタ法により非線形要素を有する系の解析を行う非線形要素を含む系の解析プログラムを記録した記録媒体において、
動作点を前記非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更した特性を用いて電子計算機に計算を行なわせる特性補償法のプログラムを記憶していることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体。
【請求項６】
請求項５に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法において、変数変換により、各段における移動距離をＬ倍に拡大又は縮小する操作に対し、前記非線形要素の特性曲線を１／Ｌとすることにより動作点が非線形要素の特性曲線上に乗るように前記非線形要素の特性を変更したことを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体。
【請求項７】
請求項５または請求項６に記載する非線形要素を含む系の解析プログラムを記録した記録媒体において、
前記陰的ルンゲクッタ法は二段対角型陰的ルンゲクッタ法であることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体。
【請求項８】
請求項５〜請求項７の何れか一つに記載する非線形要素を含む系の解析プログラムを記録した記録媒体において、
前記系は電気回路であり、前記非線形要素は前記電気回路を構成する非線形素子であることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は非線形要素を含む系の解析装置および解析プログラムを記録した記録媒体に関し、特にダイオード等の強い非線形性を有する非線形素子を含む回路からなるパワーエレクトロニクス機器のシミュレーション解析等に適用して有用なものである。
【背景技術】
【０００２】
非線形要素を含む系の代表例として非線形素子を含む回路からなるパワーエレクトロニクス機器がある。この種のパワーエレクトロニクス機器の系統への導入においてはシミュレーション解析が欠かせないものとなっており、そのためのプログラムとして、優れた解析性能を持つ瞬時値解析プログラム（以下、ＸＴＡＰ（expandable Transient Analysis Program）と称する）が提案されている（非特許文献１参照）。ＸＴＡＰでは回路計算を行う際の積分手法に二段対角型陰的ルンゲクッタ法（以下、２Ｓ−ＤＩＲＫともいう）を用いている。２Ｓ―ＤＩＲＫとは、対角型陰的ルンゲクッタ法として分類される数値積分手法の一種で、図１６に示すように、Ａ）１段目の計算、Ｂ）変数変換、Ｃ）２段目の計算、の３つの計算結果を経て、時間刻みｈ後の値を計算する。時間刻みｈ進むために、２回の積分計算と１回の変数変換を行う。ここで、変数変換（２段目の開始点を求めるための変換）とは、開始点から点Ａまでの距離を１＋√２倍、線形に拡大する変換をいう。
【０００３】
ところが、２Ｓ−ＤＩＲＫでダイオードを含む回路など非線形性の強い回路を解析すると、電圧や電流に数値計算上のスパイクが発生する場合がある。このスパイクは現実には発生しないため、使用者が実際の過電圧、過電流と区別する必要がある。
【０００４】
かかるスパイクの発生という問題は台形則でも発生する事象であり、非特許文献２をはじめ、いくつかの対処法が既に提案されている。具体的には、１）スナバを非線形素子に接続する（非特許文献３参照）、２）後退オイラー法を用いる、３）電流零クロス点を見つけて補間する（零点補間）等である。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】電力中央研究所 研究報告Ｈ０６００２
【非特許文献２】パワーエレクトロニクスシステムにおけるモデリングとシミュレーション技術協同研究委員会 編，「パワーエレクトロニクスシステムにおけるモデリングとシミュレーション技術」，電気学会技術報告１１１４号、（２００８）
【非特許文献３】P. kuffel, K. Kent and G. Irwin, “The Implementation and Effectiveness of Linear Interpolation within Digital Simulation,”International Conf. on Power Systems Transients IPST pp. ４９９−５０４、（１９９５）
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、スナバを非線形素子に接続する場合、スナバの値を適切に設計する必要があることや、スナバを接続することで計算負荷が増大することから、シミュレーションの手間が増大する。それ以外の方法も、スイッチングを正確に検出しなければいけないことや、数値振動が発生してしまった後は対処が難しい等、新たな課題も多い。
【０００７】
本発明は、上記問題点に鑑み、陰的ルンゲクッタ法、特に２Ｓ−ＤＩＲＫや３Ｓ−ＤＩＲＫを含む対角型陰的ルンゲクッタ法として分類される数値積分手法を利用してシミュレーション解析を行う際のノイズとなるスパイクの発生を防止し得る非線形要素を含む系の解析装置および解析プログラムを記録した記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成する本発明は、次の知見を基礎とするものである。
図１は直流電圧源とリアクトルとダイオードが直列接続された回路であり、ここでは例題回路と呼ぶ。リアクトルには初期電流ｉ０としてダイオード順方向の電流が与えられており、ダイオードはＯＮとなっている。電源電圧が、電流方向とは逆方向のため、一定時間後、リアクトルに流れる電流は０となりダイオードはＯＦＦとなるが、そのとき初期電流の値によっては２Ｓ−ＤＩＲＫで計算したダイオードやリアクトルの電圧にスパイクが生じる場合がある。
【０００９】
このスパイクは解析的には発生しない現象であり、あくまで数値計算上発生するものである。このとき、スパイクが発生するかどうかは初期電流ｉ０に依存する。また、スパイクが発生する場合でもその大きさは、同じく初期電流ｉ０に依存する。
【００１０】
図２に様々な積分手法で例題回路をシミュレートした場合の計算結果を示す。同図（ａ）に示すような理想解の特性に対し、同図（ｂ）に示す２Ｓ−ＤＩＲＫでは、数値計算的なスパイクが発生している。同図（ｃ）に示す台形則ではスパイクが発生する上に、それが持続する。一方、図２（ｂ）に示す２Ｓ−ＤＩＲＫでは、台形則のように持続振動は発生しない。
【００１１】
２Ｓ−ＤＩＲＫで、スパイクが発生する理由を説明するため、ダイオードがＯＮからＯＦＦに変化する際の模式図を図３に示す。現象を分かりやすくするため、図３では１段目終了時の値や後述する式（３）による変数変換後の状態も図示したが、通常シミュレーション結果として採用されるのは、２段目終了時（図４の点Ｃ）の値のみである。
【００１２】
非線形素子の収束計算は１段目と２段目の終了時でしか行われないため、変数変換後には電流が負になっているのにもかかわらず、ダイオードはＯＮと判定されている。このことに由来する行き過ぎが、変数変換後（つまり２段目計算開始時）に発生しており、２段目の計算は、動作点が素子の特性から外れた点から開始することとなる。そして、２段目開始時においてはダイオードに負の電流が流れているため、ダイオードはＯＦＦと判定され、ダイオードに流れる電流は０となる。このとき、リアクトルの電流を負から急速に０にするためスパイク電圧が発生している。
【００１３】
ここで、現象を整理するため、
Ｃａｓｅ１.  １段目で負電流となりダイオードがＯＦＦする場合、
Ｃａｓｅ２.  １段目では正電流だが、変数変換後に負電流となる場合、
Ｃａｓｅ３.  １段目、変数変換後ともに正電流で、２段目に負電流となる場合、と場合分けをしたときの計算過程を示す。
【００１４】
＜Ｃａｓｅ１＞  このときの計算過程を図４に示す。ダイオードは１段目の時点でＯＦＦと判定される。この結果、ダイオードは大抵抗となるが、変数変換後の電流は負電流となる。
【００１５】
図４に示す通り、このときのダイオード電圧、電流は本来のダイオードの特性から外れており、２段目の計算でリアクトルの負電流を急速に０にするためスパイク電圧が発生する。
【００１６】
＜Ｃａｓｅ２＞  このときの計算過程を図５に示す。ダイオードは１段目の時点でＯＮと判定されているが、変数変換後の電流は負電流となっている。２段目の計算ではダイオードはＯＦＦと判定される。この結果、ダイオードは大抵抗となり、リアクトルの負電流を急速に０にするため、この場合もスパイク電圧が発生する。１段目計算完了時点から変数変換完了後の間、ダイオードはＯＮとして判定されていることとなる。このケースでも、変数変換後（２段目開始時）にダイオード電圧、電流は本来のダイオードの特性から外れている。
【００１７】
＜Ｃａｓｅ３＞  このときの計算過程を図６に示す。ダイオードは１段目の時点でＯＮと判定されており、変数変換後の電流も正電流である。ダイオードは２段目の計算中にＯＮからＯＦＦへと変化するが、この場合スパイク電圧は発生しない。これは、ダイオードの特性が自身の非線形特性から外れておらず、想定どおりの動作をしているからである。
【００１８】
このように、変数変換を行うことにより回路の動作点が素子の特性から外れてしまう場合に、スパイク電圧が発生するという知見を得た。
【００１９】
なお、このような数値計算に由来するスパイクは、ダイオードのような非線形性の強い素子で顕著に現れるが、変数変換を行うことにより、回路の動作点が素子の特性から外れてしまう問題は、２Ｓ−ＤＩＲＫを積分手法に選択するに当たって非線形素子では常に発生する。特性が非線形である場合、図７に示すように、変数変換を行うことで動作点は常に特性から外れることとなる。ＸＴＡＰでは、区分折れ線近似を用いて非線形性を表現しているため、素子が１段目または変数変換中に折れ線区分をまたぐ場合に、動作点が素子特性から外れてしまうこととなる。区分折れ線の傾きの変化が小さく、時間刻みが小さいときは、特性のずれはわずかであるが、理想ダイオードなど急激に特性の傾きが変化する場合、このずれが顕著となる。
【００２０】
変数変換により回路の動作点が、素子の特性から外れてしまう場合にスパイクが発生することが分かった。したがって、変数変換後も動作点が素子の特性上に存在するように計算過程を変更すれば、スパイク電圧は発生しないと考えられる。
【００２１】
かかる知見を基礎とする本発明の第１の態様は、
陰的ルンゲクッタ法により非線形要素を有する系の解析を行う非線形要素を含む系の解析装置であって、
動作点を前記非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更した特性を用いて計算を行う特性補償手段を有することを特徴とする非線形要素を含む系の解析装置にある。
【００２２】
本態様によれば、当該系のシミュレーション解析に伴うスパイクの発生が、陰的ルンゲクッタ法における変数変換の際の非線形要素の特性上の行き過ぎにあるという知見に鑑み、動作点を前記非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更したので、前記スパイクの発生を良好に防止することができる。
【００２３】
本発明の第２の態様は、
第１の態様に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法において、変数変換により、各段における移動距離をＬ倍に拡大又は縮小する操作に対し、前記非線形要素の特性曲線を１／Ｌとすることにより動作点が非線形要素の特性曲線上に乗るように前記非線形要素の特性を変更したことを特徴とする非線形要素を含む系の解析装置にある。
【００２４】
本態様によれば、動作点が非線形要素の特性曲線上に確実に乗るので、最も効率的にスパイクの発生等を抑制し得る。
【００２５】
本発明の第３の態様は、
第１または第２の態様に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法は二段対角型陰的ルンゲクッタ法であることを特徴とする非線形要素を含む系の解析装置にある。
【００２６】
本態様によれば、陰的ルンゲクッタ法を二段対角型陰的ルンゲクッタ法としたので、計算負荷をさほど大きくすることなく所定の数値計算を高精度に実行し得る。
【００２７】
本発明の第４の態様は、
第１〜第３の態様のいずれか一つに記載する非線形要素を含む系の解析装置において、
前記系は電気回路であり、前記非線形要素は前記電気回路を構成する非線形素子であることを特徴とする非線形要素を含む系の解析装置にある。
【００２８】
本態様によれば、ダイオードやトランジスタ等、一般的に多くの非線形素子を有する電気回路のシミュレーション解析に最適なものとなる。
【００２９】
本発明の第５の態様は、
陰的ルンゲクッタ法により非線形要素を有する系の解析を行う非線形要素を含む系の解析プログラムを記録した記録媒体において、
動作点を前記非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更した特性を用いて電子計算機に計算を行なわせる特性補償法のプログラムを記憶していることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体にある。
【００３０】
本態様によれば、動作点を非線形要素の特性曲線に近づけるように、各段の計算に関する前記非線形要素の特性を所定量変更しているので、電子計算機によるシミュレーション結果におけるスパイク等の発生を良好に防止することができる。
【００３１】
本発明の第６の態様は、
第５の態様に記載する非線形要素を含む系の解析装置において、
前記陰的ルンゲクッタ法において、変数変換により、各段における移動距離をＬ倍に拡大又は縮小する操作に対し、前記非線形要素の特性曲線を１／Ｌとすることにより動作点が非線形要素の特性曲線上に乗るように前記非線形要素の特性を変更したことを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体にある。
【００３２】
本態様によれば、動作点が非線形要素の特性曲線上に確実に乗るので、最も効率的にスパイクの発生等を抑制し得るシミュレーションを電子計算機に実行させることができる。
【００３３】
本発明の第７の態様は、
第５または第６の態様に記載する非線形要素を含む系の解析プログラムを記録した記録媒体において、
前記陰的ルンゲクッタ法は二段対角型陰的ルンゲクッタ法であることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体にある。
【００３４】
本態様によれば、陰的ルンゲクッタ法を二段対角型陰的ルンゲクッタ法としたので、シミュレーションの計算負荷をさほど大きくすることなく所定の数値計算を電子計算機が高精度に実行し得る。
【００３５】
本発明の第８の態様は、
第５〜第７の態様の何れか一つに記載する非線形要素を含む系の解析プログラムを記録した記録媒体において、
前記系は電気回路であり、前記非線形要素は前記電気回路を構成する非線形素子であることを特徴とする非線形要素を含む系の解析プログラムを記録した記録媒体にある。
【００３６】
本態様によれば、ダイオードやトランジスタ等、一般的に多くの非線形素子を有する電気回路のシミュレーション解析を電子計算機に適切に実行させることができる。
【発明の効果】
【００３７】
本発明によれば、演算時間の増大を招来することなく、高精度のシミュレーション解析を行うことが可能となる。すなわち、特性補償法を用いることで、演算時間の増大を招来することなく、非線形性の強い系であってもスパイクが発生することを防止できる。さらに具体的には、スイッチの動作判定等は必要なく、自動的に特性の改善が行われるとともに、数値振動が発生した場合も、すぐに減衰する等の特徴から簡便かつ数値安定性が高く、実用性も高い。したがって、従来の台形則等よりも簡便かつ数値安定性の高い手法を提供することができ、特に非線形性が強いパワーエレクトロニクス機器を含んだ回路の解析に有効である。
【図面の簡単な説明】
【００３８】
【図１】本発明の基礎となる知見を得るためのスパイクが発生する例題回路を示す回路図である。
【図２】図１に示す例題回路を様々な積分手法でシミュレーションした場合の計算結果を示す特性図である。
【図３】２Ｓ−ＤＩＲＫで、スパイクが発生する理由を説明するため、ダイオードがＯＮからＯＦＦに変化する際の態様を示す模式図である。
【図４】図１に示す例題回路において、１段目で負電流となりダイオードがＯＦＦする場合の計算過程を示す特性図である。
【図５】図１に示す例題回路において、１段目では正電流だが、変数変換後に負電流となる場合の計算過程を示す特性図である。
【図６】図１に示す例題回路において、１段目、変数変換後ともに正電流で、２段目に負電流となる場合の計算過程を示す特性図である。
【図７】従来の２Ｓ−ＤＩＲＫにおける変数変換に伴い発生するずれを説明するための特性図である。
【図８】従来の２Ｓ−ＤＩＲＫの積分過程を概念的に示す説明図である。
【図９】従来の２Ｓ−ＤＩＲＫのＸＴＡＰにおける計算過程を示すフローチャートである。
【図１０】本発明の実施の形態に係る解析装置を示すブロック図である。
【図１１】特性補償法の概念を示す説明図である。
【図１２】特性補償法を適用し、図１に示す例題回路において、１段目でダイオードがＯＦＦする場合の計算過程を示す特性図である。
【図１３】特性補償法を適用し、図１に示す例題回路において、１段目、変数変換後ともに正電流で、２段目でダイオードがＯＦＦする場合の計算過程を示す特性図である。
【図１４】特性補償法の適用前後の計算ケースの分類を示す説明図である。
【図１５】特性補償法を含む２Ｓ−ＤＩＲＫのＸＴＡＰにおける計算過程を示す本発明の実施の形態に係るフローチャートである。
【図１６】ＸＴＡＰにおける積分手法である２Ｓ−ＤＩＲＫの原理を説明するための説明図である。
【発明を実施するための形態】
【００３９】
以下、本発明の実施の形態を図面に基づき詳細に説明する。
２Ｓ−ＤＩＲＫにおけるスパイク発生の原因が、変数変換の行き過ぎにある点は前に述べたが、スパイク発生防止のための特性補償機能（後に詳述する）を追加した本形態における２Ｓ−ＤＩＲＫの説明とともに、本形態に係る解析装置を説明する。
【００４０】
まず、従来の２Ｓ−ＤＩＲＫを概説する。２Ｓ−ＤＩＲＫは、２次の精度を持つ２段の積分計算手法であり、一時的に中間の時刻での値を算出してから、次のステップの値を計算する。すなわち、２Ｓ−ＤＩＲＫでは、微分方程式
【数１】
2014222486.tif 000003
に対する、時間刻みｈで、時刻ｔｎ−１からｔｎに進むための数式は次のとおりとなる。
【００４１】
まず１段目の公式
【数２】
2014222486.tif 000004
を用いて、時刻ｔｎ−１からｈ~だけ進んだ時刻ｔ~ｎにおける一時的な解ｙ~ｎを算出する。
【００４２】
次に
【数３】
2014222486.tif 000005
として、２段目の公式
【数４】
2014222486.tif 000006
を用いて、時刻ｔｎにおける解ｙｎを求める。また、式（２）〜（４）は、
【数５】
2014222486.tif 000007
とも表すこともできる。
【００４３】
かかる２Ｓ−ＤＩＲＫの積分過程は図８のようにイメージできる。同図に示すように、２Ｓ−ＤＩＲＫにおいて、１段目の計算はｈ~だけ後退オイラー法で計算するものである。次に式（３）による変数変換（以下、変数変換と称す）を行い、図８のＡ点までの移動距離を（１＋√２）倍にする線形補間した後、さらにｈ~だけ後退オイラー法で計算する。
【００４４】
かかる２Ｓ−ＤＩＲＫのＸＴＡＰにおける計算過程をフローチャートにしたものを図９に示す。なお、同図中、Ｆ；各素子の接続情報やコンダクタンスの情報が格納されたタブロー行列と呼ばれる正方行列、ｘ；各素子の電圧、電流情報を格納する列ベクトル、ｙ；回路中の電圧源、電流源の値を格納する列ベクトルである。
【００４５】
同図に示すように、従来の２Ｓ−ＤＩＲＫによるＸＴＡＰでは、図中にＳＴ１〜ＳＴ８で示す各処理により所定の計算が進められる。ここで、ＳＴ３〜ＳＴ５の処理が２Ｓ−ＤＩＲＫの１段目の処理である。また、ＳＴ５、ＳＴ８における非線形素子の収束計算には主にニュートン・ラフソン法を使用しており、その収束判定は１段目終了時点と２段目終了時点でそれぞれ行っている。
【００４６】
積分手法が、上述の如き２Ｓ−ＤＩＲＫであれば、数値計算的な振動が発生してもしだいに減衰し、持続することがないことは確認されている。しかし、従来の２Ｓ−ＤＩＲＫは、数値計算的な振動の減衰は保証しても、数値計算的な振動が発生しないことは保証していない。ちなみに、２Ｓ−ＤＩＲＫを適用しても、非線形性が強いパワーエレクトロニクス回路等では解析解とは異なる数値計算的なスパイクが生じることがある。
【００４７】
そこで、本形態における２Ｓ−ＤＩＲＫでは、１段目の素子の特性を、現在の動作点を中心に１／（１＋√２）に縮小する処理を行なう。
【００４８】
図８に示すように、従来の２Ｓ−ＤＩＲＫでは、変数変換によりｔｎ−１からｔ~ｎまでの移動距離を（１＋√２）倍拡大することになるため、特性を１／（１＋√２）に縮小すれば変数変換後にも、動作点が素子の特性上に存在することとなる（この点は、後に証明する）。すなわち、変数変換後の動作点が必ず非線形特性上に乗ることを保証することができる。かかる手法を特性補償法(Characteristics Compensation Method: CC-Method)と呼ぶ。
【００４９】
例えば図１に示すような、非線形素子（ダイオード）を有する例題回路のシミュレーション解析を行うための本形態に係る解析装置は、前記特性補償法を実行処理する機能を有するものである。図１０にそのブロック図を示す。同図に示すように、解析装置Ｉは、ＲＡＭおよびＲＯＭで構築されるデータベース１，２と、データベース１，２との間で適宜情報の授受を行いながら所定の演算処理を実行する演算処理部３とを有している。データベース１，２に対する必要な情報の入力はキーボード等の入力手段４を介して行い、演算処理部３で所定の処置をして得たシミュレーション解析結果等はプリンタ等の出力手段５を介して行われる。ここで、ＲＡＭであるデータベース１には、解析対象となる回路のパラメータ等が入力手段４を介して演算に先立ち、都度入力される。ＲＯＭであるデータベース２には、演算処理部３で実行する計算手順をあらわすプログラム等が記憶されている。ここで、演算処理部３は、データベース２に記憶されているプログラムに基づき前記特性補償法を実行する特性補償手段３Ａを有している。ここで、演算処理部３で処理する演算手順を記載したプログラムは、外部の記録媒体に記録しておくこともでき、この場合には当該記録媒体が記憶しているデータを入力手段４を介してデータベース２に読み込むようにすれば良い。
【００５０】
かかる本形態の解析装置Ｉは、特性補償法を実行する特性補償手段３Ａを有するので、すなわち本形態における２Ｓ−ＤＩＲＫでは、１段目の素子の特性を、現在の動作点を中心に１／（１＋√２）に縮小する処理を行なうので、非線形素子を含む回路であっても、そのシミュレーション解析結果からスパイクの発生を除去し得る。
【００５１】
特性補償法の概念図を図１１に示す。ここで、図１の例題回路に特性補償法を適用した場合、計算過程を以下のＣａｓｅ４、Ｃａｓｅ５の２パターンに場合分けすることができる。それぞれの計算過程を図１２および図１３に示す。
【００５２】
Ｃａｓｅ４；１段目でダイオードがＯＦＦする場合（図１２に示す場合）、
Ｃａｓｅ５；１段目、変数変換後ともに正電流で、２段目でダイオードがＯＦＦする場合（図１３に示す場合）。
【００５３】
特性補償法の効果により、前述のＣａｓｅ２のように、変数変換中に動作点が折れ線区分をまたぐパターンは存在しない。つまり、図１４に示すように、前述のＣａｓｅ１、Ｃａｓｅ２は、Ｃａｓｅ４に統合されることとなる。Ｃａｓｅ５は、前述のＣａｓｅ３と同じ動作をするケース、すなわちもともとスパイクが発生しない場合である。
【００５４】
このように、特性補償法を用いることで、非線形性がある回路において、動作点が素子の特性から外れることを防ぐことができる。しかしながら、特性補償法をそのまま適用する場合、素子の現在の動作点を中心に特性を縮小するため、動作点が動くたびに素子の非線形特性を作り直す必要が出てくる。通常、回路の動作点は時々刻々と変化するため、実際には非線形特性を時間刻みｈごとに作り直すこととなり、計算負荷が大きくなる。
【００５５】
そこで、計算負荷を軽くするため、特性自体を縮小するのではなく、１段目計算開始時から終了時までの移動距離（ｙｎ−１とｙ~ｎ−１）を（１＋√２）倍拡大することで、等価的に特性の縮小を実現する。具体的なアルゴリズムは下記の通りである。
１）  後退オイラー法で時刻ｔｎ−１＋ｈ~における値を計算する（１段目の計算）。
２）  式（３）の変数変換を行う。
３）  変数変換後の値を用いて区分折れ線での動作点を判定する。このとき、折れ線の特性は元の特性（２段目の特性と同じ）を使用する。
４）  計算が収束した後、動作点が折れ線の区分を跨いだ場合は、下式(６）を用いて特性の切片の補正を行う（以下、切片補正と呼ぶ）。
【００５６】
【数６】
2014222486.tif 000008
ここで、ｘｎ−１，ｙｎ−１は、時刻ｔｎ−１における素子の動作点を、ａ，ｂは新しい動作点領域での傾きと切片を示す。かかる切片補正により等価的に特性の縮小を実現している。
【００５７】
以上の手順により、素子の非線形特性を時刻ごとに変更することなく、特性補償法を実現できる。
【００５８】
また、特性補償法を組み込んだ本形態に係るＸＴＡＰの計算処理のフローチャートを図１５に示す。同図において、図９に示す従来の２Ｓ−ＤＩＲＫを適用したフローチャートと同一の処理には同一の符号を付している。
【００５９】
図１５に示すように、従来の２Ｓ−ＤＩＲＫでは、１段目と２段目で同じ計算を行っていたが、特性補償法を採用した本形態では、１段目の収束判定に変数変換後の値を用いる点と、１段目の計算終了時に切片補正を行う点（ＳＴ１１およびＳＴ１２に示す処理）が異なっている。また、本形態における２段目の処理は、図９の１段目の処理と同様である。
【００６０】
なお、特性補償法が動作するのは、例題回路では、ダイオードの動作区分が変化する１ステップのみであり、それ以外のほぼすべての線形区間は従来通りの計算を行う。よって、特性補償法が計算精度を悪化させることはない。
【００６１】
次に、上述の如き特性補償法では、変数変換による（１＋√２）倍の拡大後に、回路の動作点が素子特性上に乗り、変数変換の際の行きすぎが発生せず、数値計算上のスパイクの発生を防止できる点に関連し、上述の如き特性補償を行うことで、動作点が非線形素子の特性から外れないことは数学的に保証されている。この点に関し詳細に説明しておく。
【００６２】
ある非線形素子の特性を、
【数７】
2014222486.tif 000009
とする。時刻ｔｎ−１から時刻t~ｎへと進むための数値積分を行っており、時刻ｔｎ−１におけるｘ、ｙの値をｘｎ−１、ｙｎ−１とする。上式（７）を動作点(ｘｎ−１、ｙｎ−１)周りで拡大縮小することを考える。
【００６３】
式（７）を(−ｘｎ−１、−ｙｎ−１) だけ平行移動させると、
【数８】
2014222486.tif 000010
【００６４】
上式（８）を原点を中心にｋ倍だけ縮小すると（２Ｓ−ＤＩＲＫの場合、ｋ＝１＋√２）、
【数９】
2014222486.tif 000011
【００６５】
上式（９）を(ｘｎ−１、ｙｎ−１)だけ平行移動させ、復元すると、
【数１０】
2014222486.tif 000012
となる。
【００６６】
特性補償法を適用した場合、２Ｓ−ＤＩＲＫの１段目の収束計算は、式（７）の特性ではなく、式（１０）の特性で行うこととなる。
【００６７】
式（１０）の特性を用いて計算した、１段目の計算結果をｘ~ｎ、ｙ~ｎとし、所定の変数変換を行うと、
【数１１】
2014222486.tif 000013
より、変数変換後のｘ座標、ｙ座標は
【数１２】
2014222486.tif 000014
となる。ここで、
【数１３】
2014222486.tif 000015
とおくと、変数変換後の解は
【数１４】
2014222486.tif 000016
となり、式（７）と同じ元の素子特性上に乗る。すなわち、本形態における特性補償を行うことで、動作点は非線形素子の特性上に乗ることとなる。
【００６８】
なお、上記実施の形態では、２Ｓ−ＤＩＲＫにおける特性補償法を解説したが、各段で非線形特性を変更し、動作点を特性上に乗せるという特性補償法の考え方は、変数変換を伴う多段法全てに適用可能である。
【００６９】
また、上記実施の形態における２Ｓ−ＤＩＲＫでは、１段目の素子の特性を、現在の動作点を中心に１／（１＋√２）に縮小する処理を行なっているが、これに限定するものではない。拡大する処理でも構わない。
【００７０】
さらに、上記実施の形態に示すように、変数変換によりｔｎ−１からｔ~ｎまでの移動距離を（１＋√２）倍拡大するのに対し、特性を１／（１＋√２）に縮小すれば、変数変換後にも、動作点を非線形素子の特性上に載せることができるが、これに限るものではない。変数変換により、各段における移動距離をＬ倍に拡大又は縮小する操作に対し、非線形要素の特性曲線をＬの逆数である１／Ｌとすれば、動作点が非線形要素の特性曲線上に確実に乗るようにすることができる。
【００７１】
本発明は、変数変換により回路等の系の動作点が非線形要素の特性から外れてしまう場合にスパイクが発生するという新たな知見を基礎とするものであるので、動作点の非線形要素の特性に対するズレを縮小することができれば、それなりのスパイク抑制効果を得ることは可能であると考えられる。したがって、上述の如く、動作点が非線形要素の特性曲線上に乗るように構成することは必須ではない。ただ、動作点が非線形要素の特性曲線上に乗るように構成した場合が最も良好にスパイク発生を抑制し得る点に関しては論を俟たない。
【００７２】
上記実施の形態では非線形素子を有するパワーエレクトロニクス機器のシミュレーション解析に用いるプログラムとしてＸＴＡＰを用いた場合を示し、ＸＴＡＰで回路計算を行う際の積分手法として２Ｓ−ＤＩＲＫを用いた場合を説明したが、解析プログラムをＸＴＡＰに限定するものではない。同様に、積分手法を２Ｓ−ＤＩＲＫに限定するものでもない。陰的ルングククッタであれば適用可能であり、特に対角型ルンゲクッタ法として分類される２Ｓ−ＤＩＲＫや３Ｓ−ＤＩＲＫであれば好適に適用することができ、同様のスパイク発生防止という効果を得る。
【００７３】
また、非線形要素、特に強い非線形要素を含み、シミュレーション解析に積分計算を必要とする系であれば、回路に限るものではない。例えば、構造体の振動のシミュレーション解析等にも好適に適用することができ、同様にスパイク発生防止という効果を得る。
【産業上の利用可能性】
【００７４】
本発明は非線形要素を有する回路や構造体の解析を行うためのシミュレーションを行う産業分野において有効に利用することができる。
【符号の説明】
【００７５】
Ｉ          解析装置
１，２      データベース
３          演算処理部
３Ａ        特性補償手段
４          入力装置
５          出力装置
【図１】
2014222486.tif 000017
【図３】
2014222486.tif 000018
【図４】
2014222486.tif 000019
【図５】
2014222486.tif 000020
【図６】
2014222486.tif 000021
【図７】
2014222486.tif 000022
【図８】
2014222486.tif 000023
【図９】
2014222486.tif 000024
【図１０】
2014222486.tif 000025
【図１１】
2014222486.tif 000026
【図１２】
2014222486.tif 000027
【図１３】
2014222486.tif 000028
【図１４】
2014222486.tif 000029
【図１５】
2014222486.tif 000030
【図１６】
2014222486.tif 000031
【図２】
2014222486.tif 000032