非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置
【課題】非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる方法等を提供する。
【解決手段】解ｘ（ｎ）が属する非線形要素の動作区間を特定する動作区間特定処理と、動作区間における非線形要素の直線の傾き及び切片を計算し、傾き及び切片から非線形システムを表す方程式を線形化した式に基づいて、解ｘ(ｎ)を求める処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とする処理とを備え、動作区間特定処理では、Ｘ軸を基準としたときの解ｘ（ｎ）が属する第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する動作区間に近い方を非線形要素の動作区間として特定する。
【解決手段】解ｘ（ｎ）が属する非線形要素の動作区間を特定する動作区間特定処理と、動作区間における非線形要素の直線の傾き及び切片を計算し、傾き及び切片から非線形システムを表す方程式を線形化した式に基づいて、解ｘ(ｎ)を求める処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とする処理とを備え、動作区間特定処理では、Ｘ軸を基準としたときの解ｘ（ｎ）が属する第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する動作区間に近い方を非線形要素の動作区間として特定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置に関する。
【背景技術】
【０００２】
非線形要素の特性が区分折れ線近似されたシステムの一例として、電力システム（電力系統）を例に挙げて説明する。なお、本発明は、電力システムに限らず、電気・電子回路を始めとするあらゆるシステムに適用することができるものである。
【０００３】
近年、電力システムには、動力負荷、照明負荷、情報機器負荷など、負荷の多くがパワーエレクトロニクス（以下、ＰＥと略記する。）回路を用いて接続されている。また、自然エネルギー発電設備や電力貯蔵装置もＰＥ回路により電力システムに接続される場合がほとんどであり、また、ＰＥ技術を適用した系統安定化装置の導入も始まっている。
【０００４】
ＰＥ回路はスイッチングデバイスを用いて高速に電流を入切することを動作原理とする。従って、ＰＥ回路を含む系統の解析には、従来の実効値レベルの解析だけでなく、波形レベルの解析、すなわち、瞬時値解析（他分野では過渡現象解析やトランジエント解析とも呼ばれる。）が必要となる。また、ＰＥ回路に関連する解析だけでなく、近年その重要性が認識されてきた電力品質に関する解析にも瞬時値解析が主に用いられる。さらに、従来より瞬時値解析を適用してきた各種過電圧、過電流、異常共振に関する解析についても、設備リプレースに伴う設計見直しを機にこれらの解析を実施する機会が増えている。
【０００５】
瞬時値解析の一例としては、解析対象の電気回路について回路方程式を作成し、この回路方程式をNewton-Raphson法（以下、ＮＲ法とも略記する。）を用いて線形化して得られた連立一次方程式を解いて電気回路の動作点を求める方法がある（例えば、特許文献１参照。）。なお、非特許文献１には、非線形要素の特性が区分折れ線近似された電子回路にＮＲ法を適用した旨が記載されている。
【０００６】
しかしながら、ＮＲ法は、初期値によっては収束しない場合があり、特に、電気回路など、ＩＧＢＴ（Insulted Gate Bipolar Transistor）等の非線形素子を含むような非線形性の強い回路に係る回路方程式に適用する際には、この収束性の問題が顕著となる。一方、全ての非線形素子の特性が単調増加で区分折れ線近似されている場合には必ず収束するKatzenelson法（例えば、非特許文献２参照。）という手法があるが、収束に必要な反復回数が多く、実用的でないことが知られている。
【０００７】
なお、このような問題は、非線形性の強い電気・電子回路に係る回路方程式にＮＲ法を適用する場合のみではなく、一般に、非線形要素の特性が区分折れ線近似された非線形システムについても同様に存在する。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００２−１９７４０１号公報
【非特許文献】
【０００９】
【非特許文献１】L.O. Chua, “Efficient computer algorithms for piecewise-linear analysis of resistive nonlinear networks,” IEEE Trans., Circuit Theory, vol.CT-18, no.1, pp 73-85, Jan. 1971.
【非特許文献２】J.Katzenelson, “An Algorithm for solving nonlinear resistor networks,” Bell Syst. Tech. J., pp. 1605-1620, Oct. 1965.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
本発明は、上記従来技術に鑑み、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明は、非線形要素の特性が区分折れ線近似された非線形システムの動作点を反復計算により求めるが、以降の説明では、各反復計算ステップ（求解処理）での解を単に「解」と称し、反復計算（収束判定処理）の結果、最終的に求めようとしている解を「求める解」と称する。
【００１２】
上記目的を達成するための本発明の第１の態様は、非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法であって、非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び前記非線形システムを表す方程式の初期解ｘ（０）を入力する初期設定処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間を特定する動作区間特定処理と、前記動作区間における前記非線形要素の直線の傾き及び切片を計算し、当該傾き及び切片から前記非線形システムを表す方程式を線形化した（ａ）式
【数１】
のＦ（ｎ−１）の値及びｙ（ｎ−１）の値を計算し、当該Ｆ（ｎ−１）の値及びｙ（ｎ−１）の値を用いて（ａ）式を解いて解ｘ(ｎ)を求める求解処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とし、一致しない場合は前記動作区間特定処理及び前記求解処理を再実行させる収束判定処理とを備え、前記動作区間特定処理では、解ｘ（ｎ）（ただし、ｎは１以上）の動作区間を特定する際には、Ｘ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する前記非線形要素の動作区間に近い方を前記非線形要素の動作区間として特定することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１３】
かかる第１の態様では、解ｘ（ｎ）が属する動作区間を特定する際に、一方の軸のみならず、Ｘ軸Ｙ軸の両方について動作区間を求め、それらのうち、解ｘ（ｎ−１）の動作区間に近い方を動作区間として特定する。すなわち、本発明は、従来のＮＲ法に比べて１回の反復で飛び越せる動作区間数が少ない。このため、本発明は、非線形要素の特性によっては求める解が属する動作区間からかけ離れた動作区間に解が移動してしまうということを回避することができ、収束性を高めることができる。
【００１４】
本発明の第２の態様は、第１の態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、前記初期設定処理実行後、Newton-Raphson法を用いて解を求め、所定の反復回数以内に解が求まらない場合に、前記動作区間特定処理から前記収束判定処理を実行することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１５】
かかる第２の態様では、第１の態様で述べたように両軸についての動作区間を求めることで収束性を高めつつ、計算時間の少ない従来ＮＲ法を併用することで全体的に計算時間を短縮することができる。
【００１６】
本発明の第３の態様は、第２の態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、前記収束判定処理では、所定の回数以内に収束しない場合に、Katzenelson法に切り替えて解を求めることを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１７】
かかる第３の態様では、収束する保証のあるKatzenelson法を併用することで確実に解が得られることを担保すると共に、Katzenelson法のみを用いた場合よりも計算時間を短縮することができる。
【００１８】
本発明の第４の態様は、第１〜第３の何れか一つの態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を電子計算機に実行させるプログラムにある。
【００１９】
本発明の第５の態様は、第４の態様に記載するプログラムを記録した電子計算機で読み取り可能な記録媒体にある。
【００２０】
本発明の第６の態様は、第１〜第３の何れか一つの態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するように構成されたコンピュータを有することを特徴とするシミュレーション装置にある。
【００２１】
かかる第４〜第６の態様では、いずれも、これらを用いて上記非線形システムの動作点を求めるに際し、計算時間を短縮するとともに、収束性を向上することができる。
【００２２】
本発明の第７の態様は、第６の態様に記載するシミュレーション装置において、前記非線形システムは、非線形素子を含む仮想的な電気・電子回路であり、前記電気・電子回路の制御を行う前記制御装置からの制御信号がＡ／Ｄコンバータを介して前記電気・電子回路のシミュレーションに用いられる入力データとして入力され、前記電気・電子回路の所定部分の電圧又は電流に関する計算値をＤ／Ａコンバータによりその電圧又は電流に変換して前記制御装置に出力し、前記制御装置の動作に合わせて実時間に同期させて前記電気・電子回路の動作を計算するように構成されたことを特徴とするシミュレーション装置にある。
【００２３】
かかる第７の態様では、電気・電子回路の動作点を計算するにあたり、少ない計算時間で、かつ、高い収束性で解を得ることができるので、実機である制御装置等の動作をリアルタイムで検証することができる。
【発明の効果】
【００２４】
本発明によれば、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置が提供される。
【図面の簡単な説明】
【００２５】
【図１】電子回路Ｉの回路図及び非線形抵抗の電圧−電流特性を示す図である。
【図２】本発明の実施形態１に係るシミュレーション装置の構成を示すブロック図である。
【図３】本発明の実施形態１に係るシミュレーション装置に係るコンピュータのハードウェア構成を示すブロック図である。
【図４】本発明の実施形態１に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法のフローを示す図である。
【図５】電子回路Ｉの電圧−電流特性を示す図である。
【図６】従来ＮＲ法における動作区間選定方法を示す図である。
【図７】従来ＮＲ法の収束性を示す図である。
【図８】本発明の実施形態１に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法の収束性を示す図である。
【図９】Katzenelson法の収束性を示す図である。
【図１０】本発明の実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法のフローを示す図である。
【図１１】電子回路ＩＩ及び電子回路ＩＩＩに係る非線形抵抗の電圧−電流特性を示す図である。
【図１２】電子回路ＩＩにおける反復計算の軌跡を示す図である。
【図１３】電子回路ＩＩＩにおける反復計算の軌跡を示す図である。
【図１４】太陽光発電設備のパワーコンディショナ回路図である。
【図１５】ダイオードとＩＧＢＴのモデルを示す図である。
【図１６】太陽光発電設備のパワーコンディショナ回路に係る各素子の非線形特性を示す図である。
【図１７】本発明の実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法により求められたインバータ出力電圧・電流波形である。
【図１８】１周期の反復回数の変化を示す図である。
【図１９】ＩＧＢＴの電圧波形と反復回数の関係を示す図である。
【図２０】本発明と従来ＮＲ法とでそれぞれ計算したインバータ出力電圧波形を比較した図である。
【図２１】本発明の実施形態３に係るシミュレーション装置の構成を示すブロック図である。
【発明を実施するための形態】
【００２６】
〈実施形態１〉
本発明の実施形態に係る「非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法（以下、動作点求解方法とも記載する）」が対象とする非線形システムは、非線形の連立方程式で表すことができる系であり、一般に、次の（１）式で導出できる。
【００２７】
【数２】
【００２８】
ｘは求めるべきＮ個の未知変数からなるベクトル、ｆ（・）はＮ本の非線形方程式を与えるベクトル関数である。（１）式は非線形の連立方程式であるから、これを解くには反復計算が必要となる。（１）式をｘについて線形化すれば、反復計算の（２）式が得られる。
【００２９】
【数３】
【００３０】
変数右肩の括弧内の数字は、反復ステップを示し、Ｆはｆ（・）をｘに対して偏微分して得られるＮ×Ｎ行列である。初期値ｘ（０）からスタートして、逐次、（２）式の連立一次方程式を解いていくことにより更新された解ｘ（ｎ）を得る。ただし、Ｆ（ｎ−１），ｙ（ｎ−１）は、１つ前の反復ステップにおける解ｘ（ｎ−１）から算出することができる。解ｘ（ｎ）が収束すれば、その値を求める解とし、反復計算を終了する。
【００３１】
実際には、（１）式を構成することはなく、直接（２）式を構成して反復計算を行うことが普通である。例えば電気・電子回路の場合、（２）式を構成する一つの方法としてスパースタブロー法がある。もちろん、（１）式を（２）式の形式で構成することができれば、スパースタブロー法以外の手法でもよい。
【００３２】
本発明に係る動作点求解方法は、非線形システムに含まれる非線形要素の特性が区分折れ線近似されていることを前提とし、上記反復計算により非線形システムを表す方程式を解く。なお、非線形要素の特性が区分折れ線近似されているとは、非線形要素の特性が複数の直線から構成されていることをいう。
【００３３】
以降、非線形要素の特性が区分折れ線近似されたシステムの一例として、電子回路を例に挙げ説明する。図１に、非線形抵抗を含む電子回路Ｉと、非線形要素の一例である非線形抵抗の特性（ｖ−ｉ特性）を示す。
【００３４】
図１（ａ）に示すように、非線形システムの一例である電子回路Ｉは、直流電源、抵抗、及び非線形要素の一例である非線形抵抗から構成されている。また、図１（ｂ）には、非線形抵抗ＲＮのｖ−ｉ特性が示されている。この非線形抵抗ＲＮのｖ−ｉ特性は、区分折れ線近似されている。すなわち、ｖ−ｉ特性が複数の直線から構成されている。各直線の電圧軸ｖ方向及び電流軸ｉ方向の範囲を動作区間という。図１（ｂ）の例では、非線形抵抗ＲＮのｖ−ｉ特性は、動作区間１（ｖ＜−１、ｉ＜−１）における直線部分と、動作区間２（−１≦ｖ＜１、−１≦ｉ＜１）における直線部分と、動作区間３（１≦ｖ、１≦ｉ）における直線部分とから構成されている。
【００３５】
このような非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法及び当該方法を実行するよう構成されたシミュレーション装置について詳細に説明する。
【００３６】
シミュレーション装置１は、図２に示すように、後述する初期設定処理、動作区間特定処理、求解処理及び収束判定処理を実行する各手段１１〜１４を備えたコンピュータ２０を有する装置である。図３に示すように、コンピュータ２０は、ＣＰＵ２１、ＲＡＭ２２、ＲＯＭ２３、ハードディスク２４等の記憶手段を備え、キーボードや記憶媒体読み取り装置などの入力装置３０と、ディスプレイ、プリンタ等の出力装置４０が接続されている。
【００３７】
図２に示すように、コンピュータ２０が実行可能な各手段１１〜１４は、後述する動作点求解方法における各処理を実行する手段であり、具体的には、初期設定手段１１と、動作区間特定手段１２と、求解手段１３と、収束判定手段１４とを有している。これらの各手段１１〜１４は、動作点求解方法をコンピュータ２０に実行させるためのプログラム１０として作成したものを読み込んで実行可能としたものである。なお、これらの各手段１１〜１４は電子回路として実装されていても良い。
【００３８】
ここで、図４を用いて、動作点求解方法の処理を説明する。図示するように、まず、非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び非線形システムを表す方程式の初期解ｘ（０）（反復子ｎ＝０）を入力する初期設定処理を行う（ステップＳ１）。非線形要素の数は、１ないし複数であり、非線形要素ごとの特性情報を入力する。また、この初期設定処理では、回路の構成に応じて電源などの他の要素についての特性情報も入力する。
【００３９】
この処理は初期設定手段１１により実行され、具体的には、ハードディスク２４に記録された特性情報及び初期解ｘ（０）、又は、入力装置３０を介して特性情報及び初期解ｘ（０）をＲＡＭ２２に記録する。特性情報は、各動作区間の電圧値ｖ及び電流値ｉの範囲と、各動作区間の直線を表す傾き及び切片からなり、動作区間が特定されれば、その動作区間の直線の傾きと切片が得られるようにＲＡＭ２２に記録されている。
【００４０】
次に、動作区間特定手段１２が、解ｘ（ｎ）（最初の動作区間特定処理の実行時では反復子ｎは０）が属する各非線形要素の動作区間を特定する動作区間特定処理を行う（ステップＳ２）。この処理により、非線形要素ごとに、非線形要素の複数の動作区間のうち、次の求解処理で用いる動作区間が特定されることになる。なお、動作区間特定処理についての詳細は後述する。
【００４１】
次に、求解手段１３が反復子ｎをインクリメントし、（２）式を解いて解ｘ（ｎ）（最初の求解処理の実行時では反復子ｎは１）を求める求解処理を行う（ステップＳ３）。詳言すると、動作区間特定処理において特定された動作区間における各非線形要素の直線の傾き及び切片を計算する。この計算は、初期設定処理で入力された特性情報から計算することができる。そして、各非線形要素の傾き及び切片から（２）式のＦ（ｎ−１）及びｙ（ｎ−１）の値を計算する。これにより得られたＦ（ｎ−１）及びｙ（ｎ−１）を用いて（２）式を解き、解ｘ（ｎ）を求め、解ｘ（ｎ）をＲＡＭ２２に記録しておく。
【００４２】
次に、収束判定手段１４が、解ｘ（ｎ）が収束したか否かを判定する収束判定処理を行う（ステップＳ４）。詳言すると、複数ある非線形要素ごとに、解ｘ（ｎ）が属する非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致するか否かを判定する。全ての非線形要素について、動作区間が一致する場合（ステップＳ５：Ｙｅｓ）、解ｘ（ｎ）を求める解とし、動作区間が一致しない非線形要素が一つでもあれば、解は収束していないとして（ステップＳ５：Ｎｏ）、動作区間特定処理（ステップＳ２）、求解処理（ステップＳ３）を再実行する。つまり、収束するまで、動作区間特定処理、求解処理、収束判定処理を繰り返し実行する。
【００４３】
以上のようにして、非線形システムを表す方程式の解を得ることができるのであるが、動作区間特定処理において、非線形要素の動作区間の特定方法により収束性や計算量が異なる。本発明に係る動作点求解方法は、この動作区間の特定方法を工夫し、従来のＮＲ法よりも高い収束性を有することとなった。以下、図１に示した電子回路Ｉを表す方程式を解く場合において、動作区間の特定方法に焦点を当てて説明するとともに、従来ＮＲ法よりも収束性が良いことを示す。
【００４４】
［従来ＮＲ方法］
図１（ａ）に示した電子回路Ｉは、スパースタブロー法、あるいは、別の定式化手法を用いた場合でも次の（３）式、（４）式に帰着できる。
【００４５】
【数４】
【００４６】
【数５】
【００４７】
ただし、ｅは電源電圧、Ｒは直列抵抗の値、ｖ，ｉは非線形抵抗の電圧、電流である。また、ｇ（・）は非線形抵抗のｖ−ｉ特性である（図１（ｂ）参照）。これらの（３）式、（４）式を重ねたものを図５に示す。（３）式、（４）式の交点が電子回路Ｉの解を与えることが分かる。
【００４８】
まず、非線形抵抗の特性情報と初期解ｘ（０）（反復子ｎ＝０）を入力し、ＲＡＭ２２に記録しておく（初期設定処理）。次に、解ｘ（ｎ）（最初の動作区間特定処理の実行時では反復子ｎは０）が属する非線形抵抗の動作区間を特定する（動作区間特定処理）。この動作区間の特定方法を、図６を用いて説明する。図６には、非線形抵抗のｖ−ｉ特性が示されており、ｖ−ｉ特性は複数の直線からなり、区分折れ線近似されている。このようなｖ−ｉ特性の電圧軸ｖに着目し、解ｘ（ｎ）が各動作区間のいずれの範囲にあるかを計算することで、解ｘ（ｎ）が属する動作区間を特定する。つまり、複数ある動作区間のうち解ｘ（ｎ）の電圧値ｖを含む動作区間を求める。この例では、解ｘ（ｎ）の電圧値ｖと、ＲＡＭ２２に記録された各動作区間とを比較することで、解ｘ（ｎ）は動作区間Ｓ１に属していることが求められる。
【００４９】
ただし、必ずしも電圧軸ｖに着目する必要はなく、電流軸ｉに着目しても良い。すなわち、複数ある動作区間のうち、解ｘ（ｎ）の電流値ｉを含む動作区間を求めても良い。いずれにせよ、従来ＮＲ法では、動作区間を特定するに際し、いずれか一方の軸のみを用いる。
【００５０】
次に、反復子ｎをインクリメントし、解ｘ（ｎ）（最初の求解処理の実行時では反復子ｎは１）を計算する（求解処理）。区分折れ線近似された非線形抵抗のｖ−ｉ特性のうち、動作区間特定処理において特定した動作区間の直線の傾き及び切片を求める。そして、この傾き及び切片から（２）式のＦ（ｎ−１）、ｙ（ｎ−１）を算出し、これらのＦ（ｎ−１）、ｙ（ｎ−１）を用いて（２）式を解き、解ｘ（ｎ）を求める。そして、解ｘ（ｎ）が属する動作区間と、解ｘ（ｎ−１）が属する動作区間とが同一であるならば収束したとして、計算を終了し、同一でないならば、動作区間特定処理から繰り返す。
【００５１】
上記の動作区間特定処理、求解処理、収束判定処理を図７を用いて視覚的に説明する。図７は、（３）式及び（４）式を表したグラフであり、（４）式は、３つの動作区間（図中、「区間１」、「区間２」、「区間３」と表記してある。）ごとの直線から構成されている。初期解ｘ（０）が予め与えられるので、解ｘ（０）が属する動作区間が定まる。図示するように、動作区間３が特定される。
【００５２】
次に、非線形抵抗のｖ−ｉ特性のうち、動作区間３の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算する。これにより、（４）式が、ａｖ＋ｂ＝ｉのように線形化され、この線形化された式と（３）式とからなる連立方程式を解いて、解ｘ（１）を得る。解ｘ（１）が属する動作区間は、動作区間１であるので、解ｘ（０）の動作区間３と異なるため、収束していないと判定し、次のステップに進む。
【００５３】
解ｘ（１）が属する動作区間を求めると、動作区間１が得られる。非線形抵抗のｖ−ｉ特性のうち、動作区間１の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算し、（４）式を線形化したａｖ＋ｂ＝ｉと（３）式とからなる連立方程式を解いて、解ｘ（２）を得る。そして、この解ｘ（２）について収束判定処理を行う。
【００５４】
以上に説明したように、解ｘ（ｎ）が求解されていくのであるが、さらに続けて、解ｘ（２）から解ｘ（３）を求めると、解ｘ（２）は初期解ｘ（０）と同じ動作区間３に属するため、解ｘ（３）は、解ｘ（１）と同じになってしまい、結局、永久に収束しない。このように、従来ＮＲ法では、解が収束しない場合がある。この非線形抵抗の場合、電流軸により動作区間を選定するようにすれば収束する。しかしながら、非線形抵抗の特性が電圧に対して電流が飽和する特性でなく（図１（ｂ）参照）、電流に対して電圧が飽和する特性の場合、やはり巡回的に解の軌跡を辿るだけで、永久に解に到達できない。
【００５５】
［動作点求解方法（両軸ＮＲ法）］
上述した従来ＮＲ法では、解ｘ（ｎ）が属する動作区間を特定するに際し、Ｘ軸（電圧軸ｖ）についてのみ考慮して動作区間を特定していたが、本発明に係る動作点求解方法では、Ｘ軸（電圧軸ｖ）及びＹ軸（電流軸ｉ）のそれぞれについて動作区間を考慮する点が異なっている。このような特徴から本発明に係る動作点求解方法を両軸ＮＲ法とも表記する。
【００５６】
図８を用いて、両軸ＮＲ法における動作点の特定方法について説明する。解ｘ（１）を求めるところまでは、従来ＮＲ法と同様である（図７参照）。解ｘ（１）から解ｘ（２）を求めるに際して、解ｘ（１）の動作区間を求めるとき、Ｘ軸（電圧軸ｖ）を基準としたときの解ｘ（ｎ）が属する非線形要素の第１動作区間を求めるとともに、Ｙ軸（電流軸ｉ）を基準としたときの解ｘ（ｎ）が属する非線形要素の第２動作区間も求める。
【００５７】
つまり、複数ある動作区間のうち、解ｘ（１）の電圧値ｖを含む動作区間を第１動作区間とし、解ｘ（１）の電流値ｉを含む動作区間を第２動作区間とする。図８に示した例では、解ｘ（１）が属する第１動作区間は区間１であり、解ｘ（１）が属する第２動作区間は区間２である。一方、解ｘ（０）が属する動作区間は、区間３であるので、この区間３に近い方である第２動作区間を動作区間として特定する。
【００５８】
以降、このように特定した動作区間（区間２）における非線形抵抗の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算する。これにより、（４）式が、ａｖ＋ｂ＝ｉのように線形化され、この線形化された式と（３）式とからなる連立方程式を解いて、解ｘ（２）を得る。この解ｘ（２）が属する動作区間は、区間２の直線上にあり、解ｘ（１）の動作区間は、先に特定したように区間２であるので、ここで収束したと判定する。
【００５９】
以上に説明したように、両軸ＮＲ法は、従来ＮＲ法では収束しないような非線形抵抗を含む電子回路に対しても収束することができる。両軸ＮＲ法は、動作区間の特定に際し、前ステップの動作区間に近い方を選択するため、従来のＮＲ法に比べて１回の反復で飛び越せる動作区間数が少ない。このため、両軸ＮＲ法は、非線形要素の特性によっては求める解が属する動作区間からかけ離れた動作区間に解が移動してしまうということを回避することができる。このような両軸ＮＲ法の特徴によれば、例えば、図１（ｂ）に示すような急峻に変化する特性を有する非線形要素を含む系に対して巡回的な軌跡に陥らずに求める解を得ることができ、従来ＮＲ法よりも収束性が良くなる。
【００６０】
［Katzenelson法］
Katzenelson法は、全ての非線形素子の特性が区分折れ線近似で表現され、かつ、単調増加の場合、いかなる初期値から出発しても必ず求める解に収束するアルゴリズムである（詳細は、非特許文献２参照）。この手法の概略を図９を用いて説明する。
【００６１】
従来ＮＲ法とKatzenelson法の差異は動作区間を選定する部分に加えて、次の反復ステップにおける解を修正する部分にある。
【００６２】
従来ＮＲ法では、現在の動作区間と新しく選定される動作区間がどれだけ離れていようとも、解が与える非線形抵抗の電圧（若しくは電流）からそのまま新しい動作区間を決定する。一方、Katzenelson法では、解が与える電圧（電流）が属する動作区間の方向に移動はするものの、最初に横切る区間と区間の境界で停止して、その向こう側を新しい動作区間とする。さらに、次の反復ステップにおける解をその境界点に修正する。すなわち、Katzenelson法は、新しい動作区間を選定する際に、一足飛びに区間を移動せずに、必ず区間の境界で停止する手法である。
【００６３】
図９に示すように、従来ＮＲ法と同様に、最初の反復ステップにおける解を求めると、点Ｃとなる。しかしながら、点Ｃを解として採用せず、初期値が属する動作区間と点Ｃ方向の隣の動作区間の境界、すなわち、区間３と区間２の境界を更新された解ｘ（１）とする。また、新たな動作区間は解ｘ（１）の向こう側に位置する区間２とする。区間２には求める解が存在するため、次の反復ステップにおいてこの解に到達し、反復計算が終了する。
【００６４】
このように、Katzenelson法は数学的に収束が保証されている。しかし、解の更新に加えて、区間の境界に位置する新たな解を計算すること、反復ステップごとに一つずつしか区間を移動できないことの理由により、従来ＮＲ法よりも計算量が増大する。また、両軸ＮＲは、従来ＮＲ法よりも計算量が増大するが、Katzenelson法のように一つずつしか区間を移動できないわけではなく、複数の動作区間を飛び越すことができるため、両軸ＮＲ法は、Katzenelson法よりも計算量は少ない。表１に、従来ＮＲ法、両軸ＮＲ法、Katzenelson法の収束性と計算量を示す。
【００６５】
【表１】
【００６６】
以上に説明した、非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法及びシミュレーション装置によれば、非線形要素の特性が区分折れ線近似された非線形システムの動作点をKatzenelson法よりも高速に、且つ従来ＮＲ法とは異なり確実に収束することができる。
【００６７】
〈実施形態２〉
非線形要素の特性が区分折れ線近似されたシステムの動作点を求めるに際し、両軸ＮＲ法のみならず、従来ＮＲ法とKatzenelson法を併用しても良い。この各方法を併用して非線形システムを表す方程式を解く方法について、図１０を用いて説明する。
【００６８】
まず、現在の反復子ｎを０で初期化する（ステップＳ１０）。次に、初期設定処理を行う（ステップＳ１１）。次に、反復子ｎをインクリメントする（ステップＳ１２）。そして、反復子ｎが閾値ＮＢ以上であるか否かを比較する（ステップＳ１３）。閾値ＮＢは、従来ＮＲ法での計算を打ち切り、両軸ＮＲ法に切り替えるステップ数を定めた所定値である。
【００６９】
反復子ｎが閾値ＮＢ未満であるならば（ステップＳ１３：Ｎｏ）、従来ＮＲ法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１４）。従来ＮＲ法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７０】
一方、反復子ｎが閾値ＮＢ以上であるときは（ステップＳ１３：Ｙｅｓ）、反復子ｎが閾値ＮＫ以上であるか否かを比較する（ステップＳ１６）。閾値ＮＫは、両軸ＮＲ法での計算を打ち切り、Katzenelson法に切り替えるステップ数を定めた所定値である。
【００７１】
反復子ｎが閾値ＮＫ未満であるならば（ステップＳ１６：Ｎｏ）、両軸ＮＲ法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１７）。両軸ＮＲ法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７２】
一方、反復子ｎが閾値ＮＫ以上であるときは（ステップＳ１６：Ｙｅｓ）、Katzenelson法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１８）。Katzenelson法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７３】
以上のように、最初は従来ＮＲ法を用いることで、短い計算時間で解を求めようとし、一定の反復回数を超えたら、従来ＮＲ法よりも高い収束性を有する両軸ＮＲ法で解を求めようとし、さらに、一定の反復回数を超えたら、Katzenelson法を用いることで、確実に解を求めようとする。このように反復回数に基づいて、従来ＮＲ法、両軸ＮＲ法、Katzenelson法を切り替えて用いることにより、非線形システムを表す方程式について、高速且つ確実に収束するように解を得ることができる。
【００７４】
〈実施例１〉
以下、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法及びシミュレーション装置を、種々の回路に適用した実施例を示す。
【００７５】
本実施例１では、二つの電子回路について説明する。２つの電子回路とも、回路構成は図１（ａ）に示したものと同じであり、非線形抵抗ＲＮの特性が異なる。図１１（ａ）に示す非線形抵抗ＲＮを含むものを電子回路ＩＩ、図１１（ｂ）に示す非線形抵抗ＲＮを含む物を電子回路ＩＩＩと表記する。
【００７６】
図１２（ａ）は、電子回路ＩＩについて従来ＮＲ法を適用したときの解ｘ（ｎ）の推移を示し、図１２（ｂ）は、電子回路ＩＩについて実施形態２に係る動作点求解方法を適用したときの解ｘ（ｎ）の推移を示す図である。
【００７７】
図１２（ａ）に示すように、従来ＮＲ法では、初期解ｘ（０）（点０）から開始した場合、点１，２，３，４，５，６，・・・を追って反復計算が行われ、解ｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）、・・・と更新されていく。点６の次は、点３に移動するため、６→３→４→５→６という無限ループに陥り、解ｘ（３）とｘ（２）の値を交互に取るだけで、永久に求める解に収束しない。
【００７８】
一方、ＮＢ−１回目までの反復では、上述の軌跡を辿るが、両軸ＮＲ法に切り替わると、１回、若しくは２回の反復で求める解に収束する。ＮＢ−１回目の反復ステップで解がｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）のいずれにあるかによって、その後の軌跡は異なるが、これら４種類の軌跡を図１２（ｂ）に示す。電子回路ＩＩでは、両軸ＮＲ法で求める解に収束したため、Katzenelson法を適用するには至らなかった。
【００７９】
図１３（ａ）は、電子回路ＩＩＩについて従来ＮＲ法を適用したときの解ｘ（ｎ）の推移を示し、図１３（ｂ）は、電子回路ＩＩについて実施形態２に係る動作点求解方法を適用したときの解ｘ（ｎ）の推移を示す図である。
【００８０】
図１３（ａ）に示すように、従来ＮＲ法では、初期解ｘ（０）（点０）から開始した場合、点１，２，３，４，５，６，７，・・・を追って反復計算が行われ、解ｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）、ｘ（４）、・・・と更新されていく。点７の次は、初期値付近に戻るため、０→１→２→・・・→６→７→０という無限ループに陥る。
【００８１】
一方、図１３（ｂ）に示すように、両軸ＮＲ法に切り替わると、１回、若しくは２回の反復で求める解に収束する。電子回路ＩＩＩにおいても、両軸ＮＲ法で求める解に収束したため、Katzenelson法を適用するには至らなかった。
【００８２】
〈実施例２〉
本実施例では、実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を、単相インバータ回路の解析に対して適用することで、実用的な解析例における当該方法の収束性能を確認する。解析に用いたインバータ回路を図１４に、また、ダイオードモデルおよびＩＧＢＴモデルの内部を図１５に示す。ダイオードは、ＳｉＣを用いたＳＢＤ（ショットキーバリアダイオード）であり、図１５に示すように静的な整流特性を非線形抵抗Ｒstで、接合容量の変化による動的な特性を非線形キャパシタンスＣjで模擬している。ＲstとＣjの非線形特性をそれぞれ図１６（ａ）、（ｂ）に示す。ＩＧＢＴはＳｉを用いたものであり、図１５（ｂ）に示すKrausのモデルで模擬している。このモデルでは、ＩＧＢＴが内包するMOSFETに対応する部分の非線形な静特性を制御電圧源Ｖmosで表現し、酸化膜が有する容量の変化を非線形キャパシタンスＣoxで模擬している。また、ベース電荷の変化により発生するベース通過電流の動特性を制御電流源Ｉqで模擬し、テール電流の動特性を制御電流源Ｉpcで模擬している。上記Ｃoxの非線形特性を図１６（ｃ）に示す。ダイオードを模擬するＤ1はオン抵抗１０-4Ω、オフ抵抗１０14Ωの理想スイッチの特性を有する非線形抵抗、Ｄ2も同様にオン抵抗１０-5Ω、オフ抵抗１０13Ωの非線形抵抗とした。なお、Ｖmos、Ｉq、Ｉpcの計算については、文献「岡田有功、菊間俊明、高崎昌洋、竹中清、小谷和也、葛巻淳彦、松本寿彰、「インバータシミュレーションプログラムの開発（その２）−実測比較による解析精度の検証−」、電力中央研究所 研究報告 R07016，2008」に詳しく記載されている。平滑用の電解コンデンサとフィルムコンデンサは、ディスクリートの部品２並列で構成されている。これらコンデンサや回路各部の配線については、実測結果から見積もった寄生インダクタンスや寄生キャパシタンスも考慮している。この回路の過渡現象をＸＴＡＰ（eXpandable Transient Analysis Programの略。電力中央研究所 研究報告 H06002，H07004，H07005，R06017，R07016を参照。）にて解析した結果は、実測結果と良好に一致する。
【００８３】
図１５では非線形特性を有する素子を点線で囲っているが、これからも分かるように本回路は多数の非線形要素を含むため、その収束が容易ではない。計算時間刻みをΔｔ＝１０ｎｓ、反復計算手法を切り替えるパラメータをＮＢ＝１００、ＮＫ＝２００に設定して、実施形態２に係る動作点求解方法により１周期分の過渡解析を行った。過渡解析であるから、計算時間ステップごと（Δｔごと）に非線形回路の解を求めるという作業を繰り返すことになる。ここで、従来ＮＲ法から両軸ＮＲ法に切り替える反復回数をＮＢ＝１００と非常に大きな値に設定しているため、両軸ＮＲ法に切り替わった場合には従来ＮＲ法が無限ループに陥って収束しなかったものと見なしてよい。
【００８４】
図１７に、実施形態２に係る動作点求解方法により求められた１周期分のインバータ出力電圧および電流の波形を示す。所期の正弦波が得られ、インバータ回路の動作が正常にシミュレーションされていることが確認できる。この計算における、各計算時間ステップでの反復回数をグラフ化したものを図１８に示す。同図より、ほとんどの計算時間ステップにおいて反復回数は６回程度であるが、所々、１００回を超えている計算時間ステップが存在することが見て取れる。反復回数が１００回を超えている計算時間ステップは、先述のように、従来ＮＲ法が収束しないケースであり、１００回目以降に反復計算手法が両軸ＮＲ法に切り替わることにより、その後数回の反復で直ちに収束している。なお、１００回という回数は、同図のグラフにおいて従来ＮＲ法から両軸ＮＲ法に切り替えたとき反復回数の差異が判別しやすくなるように敢えて大きな値を設定したものであり、実際には、１０〜２０回の反復回数で足りる。
【００８５】
どのような場合に従来ＮＲ法が収束しないのかを調べるため、左下アームおよび右下アームのＩＧＢＴ両端電圧波形と反復回数を時刻を揃えてプロットしたものを図１９に示す。同図より、従来ＮＲ法が収束しないのはＩＧＢＴ両端電圧が急変する時刻、すなわち、スイッチングの瞬間であることが分かる。ＩＧＢＴがスイッチングをしなければ、回路中の電圧・電流は急変せず、各非線形素子の動作区間も大きく動くことは無いが、スイッチングにより電圧・電流が大きく変わると各非線形素子の動作区間も大きく動くこととなり、その収束過程に無限ループが形成されれば従来ＮＲ法は収束しなくなる。本回路において、従来ＮＲ法が収束しない場合、反復計算手法が両軸ＮＲ法に切り替わった直後に少なくとも６回の反復で収束しており、前出の電子回路Ｉ〜ＩＩＩと同様、Katzenelson法が用いられることはなかった。つまり、電子回路Ｉ〜ＩＩＩおよび本インバータ回路のいずれにおいても両軸ＮＲ法の段階で収束しており、両軸ＮＲ法は非常に良い収束特性を有していると言える。また、Katzenelson法は滅多に用いられることはないが、両軸ＮＲ法が数学的に収束保証された手法でない以上、保険の意味で後段に用意されていると言える。
【００８６】
図２０に、インバータ出力電圧波形について、実施形態２に係る動作点求解方法による結果と従来ＮＲ法による結果の比較を示す。既に図１７で示したように、実施形態２に係る動作点求解方法を用いた場合、各時間ステップにおいて反復計算が正常に収束しており、問題無く１周期の計算を終了している。一方、従来ＮＲ法を用いた場合、時刻２．３ｍｓ付近で反復計算が収束せず（図２０の＋印）、解析の続行が不可能となっている。
【００８７】
以上に説明したように、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置では、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる。
【００８８】
〈実施形態３〉
上述した実施形態１及び実施形態２では、電気・電子回路について、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を、コンピュータ２０に実行させたが、このような態様に限られない。本発明の実施形態としては、例えば、コンピュータ２０で仮想的に構築した電気・電子回路と、実機とを接続し、実機の動作の検証等を行うためのリアルタイムシミュレータ装置であってもよい。
【００８９】
図２１に、本実施形態に係るシミュレーション装置の一例を示す。なお、実施形態１と同一のものには同一の符号を付し、重複する説明は省略する。
【００９０】
図示するように、コンピュータ２０は、電気・電子回路の各素子についての特性情報や初期解が予め設定されており、電気・電子回路の所定の素子や節点における電圧電流を計算するようになっている。この計算した電圧や電流の値は、Ｄ／Ａコンバータ５１により電流、電圧に変換され、電圧電流アンプ２０１で増幅されて制御装置２００に流れる。一方、制御装置２００が出力する制御信号は、Ａ／Ｄコンバータ５０により電気・電子回路における一素子の電圧、電流を表すデータとしてプログラム１０に入力される。
【００９１】
コンピュータ２０では、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するに際し、制御装置２００の動作に合わせて実時間に同期させて電気・電子回路の動作を計算する。つまり、電気・電子回路の一部をコンピュータ２０で模擬し、当該電気・電子回路に接続される制御装置２００の動作を検証することができる。
【００９２】
〈他の実施形態〉
実施形態１〜実施形態３では、いずれも非線形システムとして電気・電子回路を対象としたが、電気・電子回路に限られず、本発明は、非線形要素の特性が区分折れ線近似されたシステム全般に適用できるものである。
【００９３】
また、実施形態２では、従来ＮＲ法、両軸ＮＲ法、Katzenelson法を併用したが、必ずしも、この構成に限らず、従来ＮＲ法と両軸ＮＲ法とを併用する場合や、両軸ＮＲ法とKatzenelson法とを併用する場合など、３つの方法のうち何れか２つを併用してもよい。
【産業上の利用可能性】
【００９４】
本発明は、非線形要素の特性が区分折れ線近似されたシステムの動作点を求めることを必要とする産業分野で利用することができる。
【符号の説明】
【００９５】
１、１００ シミュレーション装置
１０ プログラム
１１ 初期設定手段
１２ 動作区間特定手段
１３ 求解手段
１４ 収束判定手段
２０ コンピュータ
２１ ＣＰＵ
２２ ＲＡＭ
２３ ＲＯＭ
２４ ハードディスク
３０ 入力装置
４０ 出力装置
【技術分野】
【０００１】
本発明は、非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置に関する。
【背景技術】
【０００２】
非線形要素の特性が区分折れ線近似されたシステムの一例として、電力システム（電力系統）を例に挙げて説明する。なお、本発明は、電力システムに限らず、電気・電子回路を始めとするあらゆるシステムに適用することができるものである。
【０００３】
近年、電力システムには、動力負荷、照明負荷、情報機器負荷など、負荷の多くがパワーエレクトロニクス（以下、ＰＥと略記する。）回路を用いて接続されている。また、自然エネルギー発電設備や電力貯蔵装置もＰＥ回路により電力システムに接続される場合がほとんどであり、また、ＰＥ技術を適用した系統安定化装置の導入も始まっている。
【０００４】
ＰＥ回路はスイッチングデバイスを用いて高速に電流を入切することを動作原理とする。従って、ＰＥ回路を含む系統の解析には、従来の実効値レベルの解析だけでなく、波形レベルの解析、すなわち、瞬時値解析（他分野では過渡現象解析やトランジエント解析とも呼ばれる。）が必要となる。また、ＰＥ回路に関連する解析だけでなく、近年その重要性が認識されてきた電力品質に関する解析にも瞬時値解析が主に用いられる。さらに、従来より瞬時値解析を適用してきた各種過電圧、過電流、異常共振に関する解析についても、設備リプレースに伴う設計見直しを機にこれらの解析を実施する機会が増えている。
【０００５】
瞬時値解析の一例としては、解析対象の電気回路について回路方程式を作成し、この回路方程式をNewton-Raphson法（以下、ＮＲ法とも略記する。）を用いて線形化して得られた連立一次方程式を解いて電気回路の動作点を求める方法がある（例えば、特許文献１参照。）。なお、非特許文献１には、非線形要素の特性が区分折れ線近似された電子回路にＮＲ法を適用した旨が記載されている。
【０００６】
しかしながら、ＮＲ法は、初期値によっては収束しない場合があり、特に、電気回路など、ＩＧＢＴ（Insulted Gate Bipolar Transistor）等の非線形素子を含むような非線形性の強い回路に係る回路方程式に適用する際には、この収束性の問題が顕著となる。一方、全ての非線形素子の特性が単調増加で区分折れ線近似されている場合には必ず収束するKatzenelson法（例えば、非特許文献２参照。）という手法があるが、収束に必要な反復回数が多く、実用的でないことが知られている。
【０００７】
なお、このような問題は、非線形性の強い電気・電子回路に係る回路方程式にＮＲ法を適用する場合のみではなく、一般に、非線形要素の特性が区分折れ線近似された非線形システムについても同様に存在する。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００２−１９７４０１号公報
【非特許文献】
【０００９】
【非特許文献１】L.O. Chua, “Efficient computer algorithms for piecewise-linear analysis of resistive nonlinear networks,” IEEE Trans., Circuit Theory, vol.CT-18, no.1, pp 73-85, Jan. 1971.
【非特許文献２】J.Katzenelson, “An Algorithm for solving nonlinear resistor networks,” Bell Syst. Tech. J., pp. 1605-1620, Oct. 1965.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
本発明は、上記従来技術に鑑み、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明は、非線形要素の特性が区分折れ線近似された非線形システムの動作点を反復計算により求めるが、以降の説明では、各反復計算ステップ（求解処理）での解を単に「解」と称し、反復計算（収束判定処理）の結果、最終的に求めようとしている解を「求める解」と称する。
【００１２】
上記目的を達成するための本発明の第１の態様は、非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法であって、非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び前記非線形システムを表す方程式の初期解ｘ（０）を入力する初期設定処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間を特定する動作区間特定処理と、前記動作区間における前記非線形要素の直線の傾き及び切片を計算し、当該傾き及び切片から前記非線形システムを表す方程式を線形化した（ａ）式
【数１】
のＦ（ｎ−１）の値及びｙ（ｎ−１）の値を計算し、当該Ｆ（ｎ−１）の値及びｙ（ｎ−１）の値を用いて（ａ）式を解いて解ｘ(ｎ)を求める求解処理と、解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とし、一致しない場合は前記動作区間特定処理及び前記求解処理を再実行させる収束判定処理とを備え、前記動作区間特定処理では、解ｘ（ｎ）（ただし、ｎは１以上）の動作区間を特定する際には、Ｘ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する前記非線形要素の動作区間に近い方を前記非線形要素の動作区間として特定することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１３】
かかる第１の態様では、解ｘ（ｎ）が属する動作区間を特定する際に、一方の軸のみならず、Ｘ軸Ｙ軸の両方について動作区間を求め、それらのうち、解ｘ（ｎ−１）の動作区間に近い方を動作区間として特定する。すなわち、本発明は、従来のＮＲ法に比べて１回の反復で飛び越せる動作区間数が少ない。このため、本発明は、非線形要素の特性によっては求める解が属する動作区間からかけ離れた動作区間に解が移動してしまうということを回避することができ、収束性を高めることができる。
【００１４】
本発明の第２の態様は、第１の態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、前記初期設定処理実行後、Newton-Raphson法を用いて解を求め、所定の反復回数以内に解が求まらない場合に、前記動作区間特定処理から前記収束判定処理を実行することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１５】
かかる第２の態様では、第１の態様で述べたように両軸についての動作区間を求めることで収束性を高めつつ、計算時間の少ない従来ＮＲ法を併用することで全体的に計算時間を短縮することができる。
【００１６】
本発明の第３の態様は、第２の態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、前記収束判定処理では、所定の回数以内に収束しない場合に、Katzenelson法に切り替えて解を求めることを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法にある。
【００１７】
かかる第３の態様では、収束する保証のあるKatzenelson法を併用することで確実に解が得られることを担保すると共に、Katzenelson法のみを用いた場合よりも計算時間を短縮することができる。
【００１８】
本発明の第４の態様は、第１〜第３の何れか一つの態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を電子計算機に実行させるプログラムにある。
【００１９】
本発明の第５の態様は、第４の態様に記載するプログラムを記録した電子計算機で読み取り可能な記録媒体にある。
【００２０】
本発明の第６の態様は、第１〜第３の何れか一つの態様に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するように構成されたコンピュータを有することを特徴とするシミュレーション装置にある。
【００２１】
かかる第４〜第６の態様では、いずれも、これらを用いて上記非線形システムの動作点を求めるに際し、計算時間を短縮するとともに、収束性を向上することができる。
【００２２】
本発明の第７の態様は、第６の態様に記載するシミュレーション装置において、前記非線形システムは、非線形素子を含む仮想的な電気・電子回路であり、前記電気・電子回路の制御を行う前記制御装置からの制御信号がＡ／Ｄコンバータを介して前記電気・電子回路のシミュレーションに用いられる入力データとして入力され、前記電気・電子回路の所定部分の電圧又は電流に関する計算値をＤ／Ａコンバータによりその電圧又は電流に変換して前記制御装置に出力し、前記制御装置の動作に合わせて実時間に同期させて前記電気・電子回路の動作を計算するように構成されたことを特徴とするシミュレーション装置にある。
【００２３】
かかる第７の態様では、電気・電子回路の動作点を計算するにあたり、少ない計算時間で、かつ、高い収束性で解を得ることができるので、実機である制御装置等の動作をリアルタイムで検証することができる。
【発明の効果】
【００２４】
本発明によれば、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置が提供される。
【図面の簡単な説明】
【００２５】
【図１】電子回路Ｉの回路図及び非線形抵抗の電圧−電流特性を示す図である。
【図２】本発明の実施形態１に係るシミュレーション装置の構成を示すブロック図である。
【図３】本発明の実施形態１に係るシミュレーション装置に係るコンピュータのハードウェア構成を示すブロック図である。
【図４】本発明の実施形態１に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法のフローを示す図である。
【図５】電子回路Ｉの電圧−電流特性を示す図である。
【図６】従来ＮＲ法における動作区間選定方法を示す図である。
【図７】従来ＮＲ法の収束性を示す図である。
【図８】本発明の実施形態１に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法の収束性を示す図である。
【図９】Katzenelson法の収束性を示す図である。
【図１０】本発明の実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法のフローを示す図である。
【図１１】電子回路ＩＩ及び電子回路ＩＩＩに係る非線形抵抗の電圧−電流特性を示す図である。
【図１２】電子回路ＩＩにおける反復計算の軌跡を示す図である。
【図１３】電子回路ＩＩＩにおける反復計算の軌跡を示す図である。
【図１４】太陽光発電設備のパワーコンディショナ回路図である。
【図１５】ダイオードとＩＧＢＴのモデルを示す図である。
【図１６】太陽光発電設備のパワーコンディショナ回路に係る各素子の非線形特性を示す図である。
【図１７】本発明の実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法により求められたインバータ出力電圧・電流波形である。
【図１８】１周期の反復回数の変化を示す図である。
【図１９】ＩＧＢＴの電圧波形と反復回数の関係を示す図である。
【図２０】本発明と従来ＮＲ法とでそれぞれ計算したインバータ出力電圧波形を比較した図である。
【図２１】本発明の実施形態３に係るシミュレーション装置の構成を示すブロック図である。
【発明を実施するための形態】
【００２６】
〈実施形態１〉
本発明の実施形態に係る「非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法（以下、動作点求解方法とも記載する）」が対象とする非線形システムは、非線形の連立方程式で表すことができる系であり、一般に、次の（１）式で導出できる。
【００２７】
【数２】
【００２８】
ｘは求めるべきＮ個の未知変数からなるベクトル、ｆ（・）はＮ本の非線形方程式を与えるベクトル関数である。（１）式は非線形の連立方程式であるから、これを解くには反復計算が必要となる。（１）式をｘについて線形化すれば、反復計算の（２）式が得られる。
【００２９】
【数３】
【００３０】
変数右肩の括弧内の数字は、反復ステップを示し、Ｆはｆ（・）をｘに対して偏微分して得られるＮ×Ｎ行列である。初期値ｘ（０）からスタートして、逐次、（２）式の連立一次方程式を解いていくことにより更新された解ｘ（ｎ）を得る。ただし、Ｆ（ｎ−１），ｙ（ｎ−１）は、１つ前の反復ステップにおける解ｘ（ｎ−１）から算出することができる。解ｘ（ｎ）が収束すれば、その値を求める解とし、反復計算を終了する。
【００３１】
実際には、（１）式を構成することはなく、直接（２）式を構成して反復計算を行うことが普通である。例えば電気・電子回路の場合、（２）式を構成する一つの方法としてスパースタブロー法がある。もちろん、（１）式を（２）式の形式で構成することができれば、スパースタブロー法以外の手法でもよい。
【００３２】
本発明に係る動作点求解方法は、非線形システムに含まれる非線形要素の特性が区分折れ線近似されていることを前提とし、上記反復計算により非線形システムを表す方程式を解く。なお、非線形要素の特性が区分折れ線近似されているとは、非線形要素の特性が複数の直線から構成されていることをいう。
【００３３】
以降、非線形要素の特性が区分折れ線近似されたシステムの一例として、電子回路を例に挙げ説明する。図１に、非線形抵抗を含む電子回路Ｉと、非線形要素の一例である非線形抵抗の特性（ｖ−ｉ特性）を示す。
【００３４】
図１（ａ）に示すように、非線形システムの一例である電子回路Ｉは、直流電源、抵抗、及び非線形要素の一例である非線形抵抗から構成されている。また、図１（ｂ）には、非線形抵抗ＲＮのｖ−ｉ特性が示されている。この非線形抵抗ＲＮのｖ−ｉ特性は、区分折れ線近似されている。すなわち、ｖ−ｉ特性が複数の直線から構成されている。各直線の電圧軸ｖ方向及び電流軸ｉ方向の範囲を動作区間という。図１（ｂ）の例では、非線形抵抗ＲＮのｖ−ｉ特性は、動作区間１（ｖ＜−１、ｉ＜−１）における直線部分と、動作区間２（−１≦ｖ＜１、−１≦ｉ＜１）における直線部分と、動作区間３（１≦ｖ、１≦ｉ）における直線部分とから構成されている。
【００３５】
このような非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法及び当該方法を実行するよう構成されたシミュレーション装置について詳細に説明する。
【００３６】
シミュレーション装置１は、図２に示すように、後述する初期設定処理、動作区間特定処理、求解処理及び収束判定処理を実行する各手段１１〜１４を備えたコンピュータ２０を有する装置である。図３に示すように、コンピュータ２０は、ＣＰＵ２１、ＲＡＭ２２、ＲＯＭ２３、ハードディスク２４等の記憶手段を備え、キーボードや記憶媒体読み取り装置などの入力装置３０と、ディスプレイ、プリンタ等の出力装置４０が接続されている。
【００３７】
図２に示すように、コンピュータ２０が実行可能な各手段１１〜１４は、後述する動作点求解方法における各処理を実行する手段であり、具体的には、初期設定手段１１と、動作区間特定手段１２と、求解手段１３と、収束判定手段１４とを有している。これらの各手段１１〜１４は、動作点求解方法をコンピュータ２０に実行させるためのプログラム１０として作成したものを読み込んで実行可能としたものである。なお、これらの各手段１１〜１４は電子回路として実装されていても良い。
【００３８】
ここで、図４を用いて、動作点求解方法の処理を説明する。図示するように、まず、非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び非線形システムを表す方程式の初期解ｘ（０）（反復子ｎ＝０）を入力する初期設定処理を行う（ステップＳ１）。非線形要素の数は、１ないし複数であり、非線形要素ごとの特性情報を入力する。また、この初期設定処理では、回路の構成に応じて電源などの他の要素についての特性情報も入力する。
【００３９】
この処理は初期設定手段１１により実行され、具体的には、ハードディスク２４に記録された特性情報及び初期解ｘ（０）、又は、入力装置３０を介して特性情報及び初期解ｘ（０）をＲＡＭ２２に記録する。特性情報は、各動作区間の電圧値ｖ及び電流値ｉの範囲と、各動作区間の直線を表す傾き及び切片からなり、動作区間が特定されれば、その動作区間の直線の傾きと切片が得られるようにＲＡＭ２２に記録されている。
【００４０】
次に、動作区間特定手段１２が、解ｘ（ｎ）（最初の動作区間特定処理の実行時では反復子ｎは０）が属する各非線形要素の動作区間を特定する動作区間特定処理を行う（ステップＳ２）。この処理により、非線形要素ごとに、非線形要素の複数の動作区間のうち、次の求解処理で用いる動作区間が特定されることになる。なお、動作区間特定処理についての詳細は後述する。
【００４１】
次に、求解手段１３が反復子ｎをインクリメントし、（２）式を解いて解ｘ（ｎ）（最初の求解処理の実行時では反復子ｎは１）を求める求解処理を行う（ステップＳ３）。詳言すると、動作区間特定処理において特定された動作区間における各非線形要素の直線の傾き及び切片を計算する。この計算は、初期設定処理で入力された特性情報から計算することができる。そして、各非線形要素の傾き及び切片から（２）式のＦ（ｎ−１）及びｙ（ｎ−１）の値を計算する。これにより得られたＦ（ｎ−１）及びｙ（ｎ−１）を用いて（２）式を解き、解ｘ（ｎ）を求め、解ｘ（ｎ）をＲＡＭ２２に記録しておく。
【００４２】
次に、収束判定手段１４が、解ｘ（ｎ）が収束したか否かを判定する収束判定処理を行う（ステップＳ４）。詳言すると、複数ある非線形要素ごとに、解ｘ（ｎ）が属する非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致するか否かを判定する。全ての非線形要素について、動作区間が一致する場合（ステップＳ５：Ｙｅｓ）、解ｘ（ｎ）を求める解とし、動作区間が一致しない非線形要素が一つでもあれば、解は収束していないとして（ステップＳ５：Ｎｏ）、動作区間特定処理（ステップＳ２）、求解処理（ステップＳ３）を再実行する。つまり、収束するまで、動作区間特定処理、求解処理、収束判定処理を繰り返し実行する。
【００４３】
以上のようにして、非線形システムを表す方程式の解を得ることができるのであるが、動作区間特定処理において、非線形要素の動作区間の特定方法により収束性や計算量が異なる。本発明に係る動作点求解方法は、この動作区間の特定方法を工夫し、従来のＮＲ法よりも高い収束性を有することとなった。以下、図１に示した電子回路Ｉを表す方程式を解く場合において、動作区間の特定方法に焦点を当てて説明するとともに、従来ＮＲ法よりも収束性が良いことを示す。
【００４４】
［従来ＮＲ方法］
図１（ａ）に示した電子回路Ｉは、スパースタブロー法、あるいは、別の定式化手法を用いた場合でも次の（３）式、（４）式に帰着できる。
【００４５】
【数４】
【００４６】
【数５】
【００４７】
ただし、ｅは電源電圧、Ｒは直列抵抗の値、ｖ，ｉは非線形抵抗の電圧、電流である。また、ｇ（・）は非線形抵抗のｖ−ｉ特性である（図１（ｂ）参照）。これらの（３）式、（４）式を重ねたものを図５に示す。（３）式、（４）式の交点が電子回路Ｉの解を与えることが分かる。
【００４８】
まず、非線形抵抗の特性情報と初期解ｘ（０）（反復子ｎ＝０）を入力し、ＲＡＭ２２に記録しておく（初期設定処理）。次に、解ｘ（ｎ）（最初の動作区間特定処理の実行時では反復子ｎは０）が属する非線形抵抗の動作区間を特定する（動作区間特定処理）。この動作区間の特定方法を、図６を用いて説明する。図６には、非線形抵抗のｖ−ｉ特性が示されており、ｖ−ｉ特性は複数の直線からなり、区分折れ線近似されている。このようなｖ−ｉ特性の電圧軸ｖに着目し、解ｘ（ｎ）が各動作区間のいずれの範囲にあるかを計算することで、解ｘ（ｎ）が属する動作区間を特定する。つまり、複数ある動作区間のうち解ｘ（ｎ）の電圧値ｖを含む動作区間を求める。この例では、解ｘ（ｎ）の電圧値ｖと、ＲＡＭ２２に記録された各動作区間とを比較することで、解ｘ（ｎ）は動作区間Ｓ１に属していることが求められる。
【００４９】
ただし、必ずしも電圧軸ｖに着目する必要はなく、電流軸ｉに着目しても良い。すなわち、複数ある動作区間のうち、解ｘ（ｎ）の電流値ｉを含む動作区間を求めても良い。いずれにせよ、従来ＮＲ法では、動作区間を特定するに際し、いずれか一方の軸のみを用いる。
【００５０】
次に、反復子ｎをインクリメントし、解ｘ（ｎ）（最初の求解処理の実行時では反復子ｎは１）を計算する（求解処理）。区分折れ線近似された非線形抵抗のｖ−ｉ特性のうち、動作区間特定処理において特定した動作区間の直線の傾き及び切片を求める。そして、この傾き及び切片から（２）式のＦ（ｎ−１）、ｙ（ｎ−１）を算出し、これらのＦ（ｎ−１）、ｙ（ｎ−１）を用いて（２）式を解き、解ｘ（ｎ）を求める。そして、解ｘ（ｎ）が属する動作区間と、解ｘ（ｎ−１）が属する動作区間とが同一であるならば収束したとして、計算を終了し、同一でないならば、動作区間特定処理から繰り返す。
【００５１】
上記の動作区間特定処理、求解処理、収束判定処理を図７を用いて視覚的に説明する。図７は、（３）式及び（４）式を表したグラフであり、（４）式は、３つの動作区間（図中、「区間１」、「区間２」、「区間３」と表記してある。）ごとの直線から構成されている。初期解ｘ（０）が予め与えられるので、解ｘ（０）が属する動作区間が定まる。図示するように、動作区間３が特定される。
【００５２】
次に、非線形抵抗のｖ−ｉ特性のうち、動作区間３の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算する。これにより、（４）式が、ａｖ＋ｂ＝ｉのように線形化され、この線形化された式と（３）式とからなる連立方程式を解いて、解ｘ（１）を得る。解ｘ（１）が属する動作区間は、動作区間１であるので、解ｘ（０）の動作区間３と異なるため、収束していないと判定し、次のステップに進む。
【００５３】
解ｘ（１）が属する動作区間を求めると、動作区間１が得られる。非線形抵抗のｖ−ｉ特性のうち、動作区間１の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算し、（４）式を線形化したａｖ＋ｂ＝ｉと（３）式とからなる連立方程式を解いて、解ｘ（２）を得る。そして、この解ｘ（２）について収束判定処理を行う。
【００５４】
以上に説明したように、解ｘ（ｎ）が求解されていくのであるが、さらに続けて、解ｘ（２）から解ｘ（３）を求めると、解ｘ（２）は初期解ｘ（０）と同じ動作区間３に属するため、解ｘ（３）は、解ｘ（１）と同じになってしまい、結局、永久に収束しない。このように、従来ＮＲ法では、解が収束しない場合がある。この非線形抵抗の場合、電流軸により動作区間を選定するようにすれば収束する。しかしながら、非線形抵抗の特性が電圧に対して電流が飽和する特性でなく（図１（ｂ）参照）、電流に対して電圧が飽和する特性の場合、やはり巡回的に解の軌跡を辿るだけで、永久に解に到達できない。
【００５５】
［動作点求解方法（両軸ＮＲ法）］
上述した従来ＮＲ法では、解ｘ（ｎ）が属する動作区間を特定するに際し、Ｘ軸（電圧軸ｖ）についてのみ考慮して動作区間を特定していたが、本発明に係る動作点求解方法では、Ｘ軸（電圧軸ｖ）及びＹ軸（電流軸ｉ）のそれぞれについて動作区間を考慮する点が異なっている。このような特徴から本発明に係る動作点求解方法を両軸ＮＲ法とも表記する。
【００５６】
図８を用いて、両軸ＮＲ法における動作点の特定方法について説明する。解ｘ（１）を求めるところまでは、従来ＮＲ法と同様である（図７参照）。解ｘ（１）から解ｘ（２）を求めるに際して、解ｘ（１）の動作区間を求めるとき、Ｘ軸（電圧軸ｖ）を基準としたときの解ｘ（ｎ）が属する非線形要素の第１動作区間を求めるとともに、Ｙ軸（電流軸ｉ）を基準としたときの解ｘ（ｎ）が属する非線形要素の第２動作区間も求める。
【００５７】
つまり、複数ある動作区間のうち、解ｘ（１）の電圧値ｖを含む動作区間を第１動作区間とし、解ｘ（１）の電流値ｉを含む動作区間を第２動作区間とする。図８に示した例では、解ｘ（１）が属する第１動作区間は区間１であり、解ｘ（１）が属する第２動作区間は区間２である。一方、解ｘ（０）が属する動作区間は、区間３であるので、この区間３に近い方である第２動作区間を動作区間として特定する。
【００５８】
以降、このように特定した動作区間（区間２）における非線形抵抗の直線の傾きと切片を求め、Ｆ（ｎ−１）、ｙ（ｎ−１）を計算する。これにより、（４）式が、ａｖ＋ｂ＝ｉのように線形化され、この線形化された式と（３）式とからなる連立方程式を解いて、解ｘ（２）を得る。この解ｘ（２）が属する動作区間は、区間２の直線上にあり、解ｘ（１）の動作区間は、先に特定したように区間２であるので、ここで収束したと判定する。
【００５９】
以上に説明したように、両軸ＮＲ法は、従来ＮＲ法では収束しないような非線形抵抗を含む電子回路に対しても収束することができる。両軸ＮＲ法は、動作区間の特定に際し、前ステップの動作区間に近い方を選択するため、従来のＮＲ法に比べて１回の反復で飛び越せる動作区間数が少ない。このため、両軸ＮＲ法は、非線形要素の特性によっては求める解が属する動作区間からかけ離れた動作区間に解が移動してしまうということを回避することができる。このような両軸ＮＲ法の特徴によれば、例えば、図１（ｂ）に示すような急峻に変化する特性を有する非線形要素を含む系に対して巡回的な軌跡に陥らずに求める解を得ることができ、従来ＮＲ法よりも収束性が良くなる。
【００６０】
［Katzenelson法］
Katzenelson法は、全ての非線形素子の特性が区分折れ線近似で表現され、かつ、単調増加の場合、いかなる初期値から出発しても必ず求める解に収束するアルゴリズムである（詳細は、非特許文献２参照）。この手法の概略を図９を用いて説明する。
【００６１】
従来ＮＲ法とKatzenelson法の差異は動作区間を選定する部分に加えて、次の反復ステップにおける解を修正する部分にある。
【００６２】
従来ＮＲ法では、現在の動作区間と新しく選定される動作区間がどれだけ離れていようとも、解が与える非線形抵抗の電圧（若しくは電流）からそのまま新しい動作区間を決定する。一方、Katzenelson法では、解が与える電圧（電流）が属する動作区間の方向に移動はするものの、最初に横切る区間と区間の境界で停止して、その向こう側を新しい動作区間とする。さらに、次の反復ステップにおける解をその境界点に修正する。すなわち、Katzenelson法は、新しい動作区間を選定する際に、一足飛びに区間を移動せずに、必ず区間の境界で停止する手法である。
【００６３】
図９に示すように、従来ＮＲ法と同様に、最初の反復ステップにおける解を求めると、点Ｃとなる。しかしながら、点Ｃを解として採用せず、初期値が属する動作区間と点Ｃ方向の隣の動作区間の境界、すなわち、区間３と区間２の境界を更新された解ｘ（１）とする。また、新たな動作区間は解ｘ（１）の向こう側に位置する区間２とする。区間２には求める解が存在するため、次の反復ステップにおいてこの解に到達し、反復計算が終了する。
【００６４】
このように、Katzenelson法は数学的に収束が保証されている。しかし、解の更新に加えて、区間の境界に位置する新たな解を計算すること、反復ステップごとに一つずつしか区間を移動できないことの理由により、従来ＮＲ法よりも計算量が増大する。また、両軸ＮＲは、従来ＮＲ法よりも計算量が増大するが、Katzenelson法のように一つずつしか区間を移動できないわけではなく、複数の動作区間を飛び越すことができるため、両軸ＮＲ法は、Katzenelson法よりも計算量は少ない。表１に、従来ＮＲ法、両軸ＮＲ法、Katzenelson法の収束性と計算量を示す。
【００６５】
【表１】
【００６６】
以上に説明した、非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法及びシミュレーション装置によれば、非線形要素の特性が区分折れ線近似された非線形システムの動作点をKatzenelson法よりも高速に、且つ従来ＮＲ法とは異なり確実に収束することができる。
【００６７】
〈実施形態２〉
非線形要素の特性が区分折れ線近似されたシステムの動作点を求めるに際し、両軸ＮＲ法のみならず、従来ＮＲ法とKatzenelson法を併用しても良い。この各方法を併用して非線形システムを表す方程式を解く方法について、図１０を用いて説明する。
【００６８】
まず、現在の反復子ｎを０で初期化する（ステップＳ１０）。次に、初期設定処理を行う（ステップＳ１１）。次に、反復子ｎをインクリメントする（ステップＳ１２）。そして、反復子ｎが閾値ＮＢ以上であるか否かを比較する（ステップＳ１３）。閾値ＮＢは、従来ＮＲ法での計算を打ち切り、両軸ＮＲ法に切り替えるステップ数を定めた所定値である。
【００６９】
反復子ｎが閾値ＮＢ未満であるならば（ステップＳ１３：Ｎｏ）、従来ＮＲ法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１４）。従来ＮＲ法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７０】
一方、反復子ｎが閾値ＮＢ以上であるときは（ステップＳ１３：Ｙｅｓ）、反復子ｎが閾値ＮＫ以上であるか否かを比較する（ステップＳ１６）。閾値ＮＫは、両軸ＮＲ法での計算を打ち切り、Katzenelson法に切り替えるステップ数を定めた所定値である。
【００７１】
反復子ｎが閾値ＮＫ未満であるならば（ステップＳ１６：Ｎｏ）、両軸ＮＲ法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１７）。両軸ＮＲ法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７２】
一方、反復子ｎが閾値ＮＫ以上であるときは（ステップＳ１６：Ｙｅｓ）、Katzenelson法を用いて、反復計算を１ステップ分（動作点特定処理〜収束判定処理）行う（ステップＳ１８）。Katzenelson法での計算結果が収束したならば（ステップＳ１５：Ｙｅｓ）、処理を終了し、収束しないならば（ステップＳ１５：Ｎｏ）、ステップＳ１２の処理に戻る。
【００７３】
以上のように、最初は従来ＮＲ法を用いることで、短い計算時間で解を求めようとし、一定の反復回数を超えたら、従来ＮＲ法よりも高い収束性を有する両軸ＮＲ法で解を求めようとし、さらに、一定の反復回数を超えたら、Katzenelson法を用いることで、確実に解を求めようとする。このように反復回数に基づいて、従来ＮＲ法、両軸ＮＲ法、Katzenelson法を切り替えて用いることにより、非線形システムを表す方程式について、高速且つ確実に収束するように解を得ることができる。
【００７４】
〈実施例１〉
以下、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法及びシミュレーション装置を、種々の回路に適用した実施例を示す。
【００７５】
本実施例１では、二つの電子回路について説明する。２つの電子回路とも、回路構成は図１（ａ）に示したものと同じであり、非線形抵抗ＲＮの特性が異なる。図１１（ａ）に示す非線形抵抗ＲＮを含むものを電子回路ＩＩ、図１１（ｂ）に示す非線形抵抗ＲＮを含む物を電子回路ＩＩＩと表記する。
【００７６】
図１２（ａ）は、電子回路ＩＩについて従来ＮＲ法を適用したときの解ｘ（ｎ）の推移を示し、図１２（ｂ）は、電子回路ＩＩについて実施形態２に係る動作点求解方法を適用したときの解ｘ（ｎ）の推移を示す図である。
【００７７】
図１２（ａ）に示すように、従来ＮＲ法では、初期解ｘ（０）（点０）から開始した場合、点１，２，３，４，５，６，・・・を追って反復計算が行われ、解ｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）、・・・と更新されていく。点６の次は、点３に移動するため、６→３→４→５→６という無限ループに陥り、解ｘ（３）とｘ（２）の値を交互に取るだけで、永久に求める解に収束しない。
【００７８】
一方、ＮＢ−１回目までの反復では、上述の軌跡を辿るが、両軸ＮＲ法に切り替わると、１回、若しくは２回の反復で求める解に収束する。ＮＢ−１回目の反復ステップで解がｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）のいずれにあるかによって、その後の軌跡は異なるが、これら４種類の軌跡を図１２（ｂ）に示す。電子回路ＩＩでは、両軸ＮＲ法で求める解に収束したため、Katzenelson法を適用するには至らなかった。
【００７９】
図１３（ａ）は、電子回路ＩＩＩについて従来ＮＲ法を適用したときの解ｘ（ｎ）の推移を示し、図１３（ｂ）は、電子回路ＩＩについて実施形態２に係る動作点求解方法を適用したときの解ｘ（ｎ）の推移を示す図である。
【００８０】
図１３（ａ）に示すように、従来ＮＲ法では、初期解ｘ（０）（点０）から開始した場合、点１，２，３，４，５，６，７，・・・を追って反復計算が行われ、解ｘ（０）、ｘ（１）、ｘ（２）、ｘ（３）、ｘ（４）、・・・と更新されていく。点７の次は、初期値付近に戻るため、０→１→２→・・・→６→７→０という無限ループに陥る。
【００８１】
一方、図１３（ｂ）に示すように、両軸ＮＲ法に切り替わると、１回、若しくは２回の反復で求める解に収束する。電子回路ＩＩＩにおいても、両軸ＮＲ法で求める解に収束したため、Katzenelson法を適用するには至らなかった。
【００８２】
〈実施例２〉
本実施例では、実施形態２に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を、単相インバータ回路の解析に対して適用することで、実用的な解析例における当該方法の収束性能を確認する。解析に用いたインバータ回路を図１４に、また、ダイオードモデルおよびＩＧＢＴモデルの内部を図１５に示す。ダイオードは、ＳｉＣを用いたＳＢＤ（ショットキーバリアダイオード）であり、図１５に示すように静的な整流特性を非線形抵抗Ｒstで、接合容量の変化による動的な特性を非線形キャパシタンスＣjで模擬している。ＲstとＣjの非線形特性をそれぞれ図１６（ａ）、（ｂ）に示す。ＩＧＢＴはＳｉを用いたものであり、図１５（ｂ）に示すKrausのモデルで模擬している。このモデルでは、ＩＧＢＴが内包するMOSFETに対応する部分の非線形な静特性を制御電圧源Ｖmosで表現し、酸化膜が有する容量の変化を非線形キャパシタンスＣoxで模擬している。また、ベース電荷の変化により発生するベース通過電流の動特性を制御電流源Ｉqで模擬し、テール電流の動特性を制御電流源Ｉpcで模擬している。上記Ｃoxの非線形特性を図１６（ｃ）に示す。ダイオードを模擬するＤ1はオン抵抗１０-4Ω、オフ抵抗１０14Ωの理想スイッチの特性を有する非線形抵抗、Ｄ2も同様にオン抵抗１０-5Ω、オフ抵抗１０13Ωの非線形抵抗とした。なお、Ｖmos、Ｉq、Ｉpcの計算については、文献「岡田有功、菊間俊明、高崎昌洋、竹中清、小谷和也、葛巻淳彦、松本寿彰、「インバータシミュレーションプログラムの開発（その２）−実測比較による解析精度の検証−」、電力中央研究所 研究報告 R07016，2008」に詳しく記載されている。平滑用の電解コンデンサとフィルムコンデンサは、ディスクリートの部品２並列で構成されている。これらコンデンサや回路各部の配線については、実測結果から見積もった寄生インダクタンスや寄生キャパシタンスも考慮している。この回路の過渡現象をＸＴＡＰ（eXpandable Transient Analysis Programの略。電力中央研究所 研究報告 H06002，H07004，H07005，R06017，R07016を参照。）にて解析した結果は、実測結果と良好に一致する。
【００８３】
図１５では非線形特性を有する素子を点線で囲っているが、これからも分かるように本回路は多数の非線形要素を含むため、その収束が容易ではない。計算時間刻みをΔｔ＝１０ｎｓ、反復計算手法を切り替えるパラメータをＮＢ＝１００、ＮＫ＝２００に設定して、実施形態２に係る動作点求解方法により１周期分の過渡解析を行った。過渡解析であるから、計算時間ステップごと（Δｔごと）に非線形回路の解を求めるという作業を繰り返すことになる。ここで、従来ＮＲ法から両軸ＮＲ法に切り替える反復回数をＮＢ＝１００と非常に大きな値に設定しているため、両軸ＮＲ法に切り替わった場合には従来ＮＲ法が無限ループに陥って収束しなかったものと見なしてよい。
【００８４】
図１７に、実施形態２に係る動作点求解方法により求められた１周期分のインバータ出力電圧および電流の波形を示す。所期の正弦波が得られ、インバータ回路の動作が正常にシミュレーションされていることが確認できる。この計算における、各計算時間ステップでの反復回数をグラフ化したものを図１８に示す。同図より、ほとんどの計算時間ステップにおいて反復回数は６回程度であるが、所々、１００回を超えている計算時間ステップが存在することが見て取れる。反復回数が１００回を超えている計算時間ステップは、先述のように、従来ＮＲ法が収束しないケースであり、１００回目以降に反復計算手法が両軸ＮＲ法に切り替わることにより、その後数回の反復で直ちに収束している。なお、１００回という回数は、同図のグラフにおいて従来ＮＲ法から両軸ＮＲ法に切り替えたとき反復回数の差異が判別しやすくなるように敢えて大きな値を設定したものであり、実際には、１０〜２０回の反復回数で足りる。
【００８５】
どのような場合に従来ＮＲ法が収束しないのかを調べるため、左下アームおよび右下アームのＩＧＢＴ両端電圧波形と反復回数を時刻を揃えてプロットしたものを図１９に示す。同図より、従来ＮＲ法が収束しないのはＩＧＢＴ両端電圧が急変する時刻、すなわち、スイッチングの瞬間であることが分かる。ＩＧＢＴがスイッチングをしなければ、回路中の電圧・電流は急変せず、各非線形素子の動作区間も大きく動くことは無いが、スイッチングにより電圧・電流が大きく変わると各非線形素子の動作区間も大きく動くこととなり、その収束過程に無限ループが形成されれば従来ＮＲ法は収束しなくなる。本回路において、従来ＮＲ法が収束しない場合、反復計算手法が両軸ＮＲ法に切り替わった直後に少なくとも６回の反復で収束しており、前出の電子回路Ｉ〜ＩＩＩと同様、Katzenelson法が用いられることはなかった。つまり、電子回路Ｉ〜ＩＩＩおよび本インバータ回路のいずれにおいても両軸ＮＲ法の段階で収束しており、両軸ＮＲ法は非常に良い収束特性を有していると言える。また、Katzenelson法は滅多に用いられることはないが、両軸ＮＲ法が数学的に収束保証された手法でない以上、保険の意味で後段に用意されていると言える。
【００８６】
図２０に、インバータ出力電圧波形について、実施形態２に係る動作点求解方法による結果と従来ＮＲ法による結果の比較を示す。既に図１７で示したように、実施形態２に係る動作点求解方法を用いた場合、各時間ステップにおいて反復計算が正常に収束しており、問題無く１周期の計算を終了している。一方、従来ＮＲ法を用いた場合、時刻２．３ｍｓ付近で反復計算が収束せず（図２０の＋印）、解析の続行が不可能となっている。
【００８７】
以上に説明したように、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法、そのプログラム、当該プログラムを記録した記録媒体及びシミュレーション装置では、非線形要素の特性が区分折れ線近似された非線形システムの動作点を高速且つ確実に収束することができる。
【００８８】
〈実施形態３〉
上述した実施形態１及び実施形態２では、電気・電子回路について、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を、コンピュータ２０に実行させたが、このような態様に限られない。本発明の実施形態としては、例えば、コンピュータ２０で仮想的に構築した電気・電子回路と、実機とを接続し、実機の動作の検証等を行うためのリアルタイムシミュレータ装置であってもよい。
【００８９】
図２１に、本実施形態に係るシミュレーション装置の一例を示す。なお、実施形態１と同一のものには同一の符号を付し、重複する説明は省略する。
【００９０】
図示するように、コンピュータ２０は、電気・電子回路の各素子についての特性情報や初期解が予め設定されており、電気・電子回路の所定の素子や節点における電圧電流を計算するようになっている。この計算した電圧や電流の値は、Ｄ／Ａコンバータ５１により電流、電圧に変換され、電圧電流アンプ２０１で増幅されて制御装置２００に流れる。一方、制御装置２００が出力する制御信号は、Ａ／Ｄコンバータ５０により電気・電子回路における一素子の電圧、電流を表すデータとしてプログラム１０に入力される。
【００９１】
コンピュータ２０では、本発明に係る非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するに際し、制御装置２００の動作に合わせて実時間に同期させて電気・電子回路の動作を計算する。つまり、電気・電子回路の一部をコンピュータ２０で模擬し、当該電気・電子回路に接続される制御装置２００の動作を検証することができる。
【００９２】
〈他の実施形態〉
実施形態１〜実施形態３では、いずれも非線形システムとして電気・電子回路を対象としたが、電気・電子回路に限られず、本発明は、非線形要素の特性が区分折れ線近似されたシステム全般に適用できるものである。
【００９３】
また、実施形態２では、従来ＮＲ法、両軸ＮＲ法、Katzenelson法を併用したが、必ずしも、この構成に限らず、従来ＮＲ法と両軸ＮＲ法とを併用する場合や、両軸ＮＲ法とKatzenelson法とを併用する場合など、３つの方法のうち何れか２つを併用してもよい。
【産業上の利用可能性】
【００９４】
本発明は、非線形要素の特性が区分折れ線近似されたシステムの動作点を求めることを必要とする産業分野で利用することができる。
【符号の説明】
【００９５】
１、１００ シミュレーション装置
１０ プログラム
１１ 初期設定手段
１２ 動作区間特定手段
１３ 求解手段
１４ 収束判定手段
２０ コンピュータ
２１ ＣＰＵ
２２ ＲＡＭ
２３ ＲＯＭ
２４ ハードディスク
３０ 入力装置
４０ 出力装置
【特許請求の範囲】
【請求項１】
非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法であって、
非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び前記非線形システムを表す方程式の初期解ｘ（０）を入力する初期設定処理と、
解ｘ（ｎ）が属する前記非線形要素の動作区間を特定する動作区間特定処理と、
前記動作区間における前記非線形要素の直線の傾き及び切片を計算し、当該傾き及び切片から前記非線形システムを表す方程式を線形化した（ａ）式
【数１】
のＦ（ｎ−１）の値及びｙ（ｎ−１）の値を計算し、当該Ｆ（ｎ−１）の値及びｙ（ｎ−１）の値を用いて（ａ）式を解いて解ｘ(ｎ)を求める求解処理と、
解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とし、一致しない場合は前記動作区間特定処理及び前記求解処理を再実行させる収束判定処理とを備え、
前記動作区間特定処理では、解ｘ（ｎ）（ただし、ｎは１以上）の動作区間を特定する際には、Ｘ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する前記非線形要素の動作区間に近い方を前記非線形要素の動作区間として特定する
ことを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項２】
請求項１に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、
前記初期設定処理実行後、Newton-Raphson法を用いて解を求め、所定の反復回数以内に解が求まらない場合に、前記動作区間特定処理から前記収束判定処理を実行することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項３】
請求項２に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、
前記収束判定処理では、所定の回数以内に収束しない場合に、Katzenelson法に切り替えて解を求める
ことを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項４】
請求項１〜請求項３の何れか一項に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を電子計算機に実行させるプログラム。
【請求項５】
請求項４に記載するプログラムを記録した電子計算機で読み取り可能な記録媒体。
【請求項６】
請求項１〜請求項３の何れか一項に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するように構成されたコンピュータを有することを特徴とするシミュレーション装置。
【請求項７】
請求項６に記載するシミュレーション装置において、
前記非線形システムは、非線形素子を含む仮想的な電気・電子回路であり、
前記電気・電子回路の制御を行う前記制御装置からの制御信号がＡ／Ｄコンバータを介して前記電気・電子回路のシミュレーションに用いられる入力データとして入力され、
前記電気・電子回路の所定部分の電圧又は電流に関する計算値をＤ／Ａコンバータによりその電圧又は電流に変換して前記制御装置に出力し、
前記制御装置の動作に合わせて実時間に同期させて前記電気・電子回路の動作を計算するように構成されたことを特徴とするシミュレーション装置。
【請求項１】
非線形要素の特性が区分折れ線近似された非線形システムの動作点を電子計算機により求める方法であって、
非線形要素の特性を複数の動作区間ごとに直線で近似した特性情報及び前記非線形システムを表す方程式の初期解ｘ（０）を入力する初期設定処理と、
解ｘ（ｎ）が属する前記非線形要素の動作区間を特定する動作区間特定処理と、
前記動作区間における前記非線形要素の直線の傾き及び切片を計算し、当該傾き及び切片から前記非線形システムを表す方程式を線形化した（ａ）式
【数１】
のＦ（ｎ−１）の値及びｙ（ｎ−１）の値を計算し、当該Ｆ（ｎ−１）の値及びｙ（ｎ−１）の値を用いて（ａ）式を解いて解ｘ(ｎ)を求める求解処理と、
解ｘ（ｎ）が属する前記非線形要素の動作区間と、解ｘ（ｎ−１）が属する動作区間とが一致した場合に解ｘ（ｎ）を求める解とし、一致しない場合は前記動作区間特定処理及び前記求解処理を再実行させる収束判定処理とを備え、
前記動作区間特定処理では、解ｘ（ｎ）（ただし、ｎは１以上）の動作区間を特定する際には、Ｘ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第１動作区間と、Ｙ軸を基準としたときの解ｘ（ｎ）が属する前記非線形要素の第２動作区間とをそれぞれ求め、第１動作区間及び第２動作区間のうち解ｘ（ｎ−１）が属する前記非線形要素の動作区間に近い方を前記非線形要素の動作区間として特定する
ことを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項２】
請求項１に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、
前記初期設定処理実行後、Newton-Raphson法を用いて解を求め、所定の反復回数以内に解が求まらない場合に、前記動作区間特定処理から前記収束判定処理を実行することを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項３】
請求項２に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法において、
前記収束判定処理では、所定の回数以内に収束しない場合に、Katzenelson法に切り替えて解を求める
ことを特徴とする非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法。
【請求項４】
請求項１〜請求項３の何れか一項に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を電子計算機に実行させるプログラム。
【請求項５】
請求項４に記載するプログラムを記録した電子計算機で読み取り可能な記録媒体。
【請求項６】
請求項１〜請求項３の何れか一項に記載する非線形要素の特性が区分折れ線近似されたシステムの動作点を電子計算機により求める方法を実行するように構成されたコンピュータを有することを特徴とするシミュレーション装置。
【請求項７】
請求項６に記載するシミュレーション装置において、
前記非線形システムは、非線形素子を含む仮想的な電気・電子回路であり、
前記電気・電子回路の制御を行う前記制御装置からの制御信号がＡ／Ｄコンバータを介して前記電気・電子回路のシミュレーションに用いられる入力データとして入力され、
前記電気・電子回路の所定部分の電圧又は電流に関する計算値をＤ／Ａコンバータによりその電圧又は電流に変換して前記制御装置に出力し、
前記制御装置の動作に合わせて実時間に同期させて前記電気・電子回路の動作を計算するように構成されたことを特徴とするシミュレーション装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【公開番号】特開２０１１−８１６７４（Ｐ２０１１−８１６７４Ａ）
【公開日】平成２３年４月２１日（２０１１．４．２１）
【国際特許分類】
【出願番号】特願２００９−２３４５４９（Ｐ２００９−２３４５４９）
【出願日】平成２１年１０月８日（２００９．１０．８）
【出願人】（０００１７３８０９）財団法人電力中央研究所 (1,040)
【Ｆターム（参考）】
【公開日】平成２３年４月２１日（２０１１．４．２１）
【国際特許分類】
【出願日】平成２１年１０月８日（２００９．１０．８）
【出願人】（０００１７３８０９）財団法人電力中央研究所 (1,040)
【Ｆターム（参考）】
[ Back to top ]