(57)【要約】  電力ネットワーク（１）のトポロジを識別するための方法は、ネットワークパラメータの情報を用いずにネットワークトポロジの自動的且つ効率的識別を可能にする。本方法は、相互電流感度係数の推定、及び推定された感度係数からネットワーク接続行列を得るためのアルゴリズムに基づく。このアルゴリズムは、特定のノードに接続された分岐における測定された電流の変動に最も感度の高い計測ユニットが、該特定のノードのすぐ上流の物理ノードに配置された計測ユニットであるという一般的仮定に基づく。本方法は、測定値におけるノイズの存在及びこのノイズの時間相関を効果的に考慮する。【選択図】図１
2020533926.tif 000002
【特許請求の範囲】
【請求項１】
電力ネットワーク（１）のトポロジを識別するための方法であって、前記電力ネットワークは、スラックノード（Ｎ１）と、他のノードのセット（Ｎ２、．．．、Ｎ３）と、前記スラックノードに接続されたスラック分岐（１０）と、前記スラックノードと前記他のノードとを互いに接続するように構成された他の分岐のセット（Ｌ１、Ｌ２）と、を備え、前記電力ネットワークは更に、計測ユニット（Ｍ１、．．．、Ｍ７）を含む監視インフラストラクチャを備えており、前記計測ユニットのうちの１つ（Ｍ１）は、前記スラック分岐（１０）を通って前記スラックノード（Ｎ１）に流入又はそれから流出する電流を測定するように構成されており、前記他の計測ユニット（Ｍ２、．．．、Ｍ７）は各々、前記他の分岐のうちの少なくとも１つを通って前記ノードのうちの１つに流入又はそれから流出する前記電流を測定するように構成され、前記監視インフラストラクチャは、通信ネットワークに接続された処理ユニット（７）を備え、前記計測ユニットは、前記処理ユニットとの間のデータ送信が可能になるように前記通信ネットワークに接続され、
前記方法は、
Ｉ．前記計測ユニット（Ｍ１、．．．、Ｍ７）に、時間ウィンドウ（τ）にわたって繰り返して電流強度値Ｉ（ｔ）及び前記電流強度値のタイムスタンプｔ∈｛ｔ1，．．．，ｔm｝を同時に測定させるステップと、
ＩＩ．前記電流強度値の各々から以前の前記電流強度値を減算することにより、前記各計測ユニットによって測定された前記電流強度の変動ΔＩi（ｔ）を計算して、前記各計測ユニット（Ｍ１、．．．、Ｍ７）によって測定された前記電流強度の変動の時系列順序付きテーブルを編集するステップと、
ＩＩＩ．実際の前記電流強度の変動（ΔＩi（ｔ））と最尤推定（ＭＬＥ）によって予測される変動との間の差異に対応する誤差項間の負の１次系列相関を考慮しながら、ステップＩＩの間に編集されたように、前記各測定ユニット（Ｍ１、．．．、Ｍ７）によって測定された電流強度の変動ΔＩi（ｔ）を、前記他の測定ユニットによって測定された前記電流強度の変動ΔＩj（ｔ）（ｊ≠ｉ）に関連付ける電流感度係数ＫＩＩi,jの前記最尤推定を実行して、前記推定された電流感度係数ＫＩＩi,jから電流感度係数行列を得るステップと、
ＩＶ．特定のノードに接続された分岐における前記測定された電流の前記変動に最も感度の高い前記計測ユニットが、前記特定のノードのすぐ上流の前記物理ノードに配置された前記計測ユニットであるという一般的仮定に基づいて、ステップＩＩＩにおいて推定された前記電流感度係数からネットワーク接続行列ＡN×(N-1)を計算するステップと、
を含む、ことを特徴とする方法。
【請求項２】
前記他の計測ユニット（Ｍ２、．．．、Ｍ７）は各々、前記他の分岐のうちの少なくとも１つを通って前記他のノードのうちの異なる１つに流入又はそれから流出する前記電流を測定するように構成される、請求項１に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項３】
前記ステップＩＶは、以下のサブステップ、すなわち、
Ａ．前記電流感度係数行列（ＫＩＩ）の各列に対して、最大絶対値を有する非対角要素（ＫＩＩi≠j）を識別し、この要素（ＫＩＩi,j）に値１を割り当て、対角要素（ＫＩＩi,j）にも値１を割り当て、更に他の全ての非対角要素に値０を割り当てるステップ、
を含む、請求項１又は２に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項４】
前記ステップＩＶは更に、
Ｂ．前記スラックノードを見つけて、前記スラック電流に対応する前記列を前記電流感度行列から削除するサブステップを含む、請求項３に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項５】
前記ステップＩＶは、以下のサブステップ、すなわち、
ａ．前記電流感度係数行列の各要素に、対応する要素の推定値の絶対値を割り当てることによって「絶対電流感度行列」を得るサブステップと、
ｂ．前記絶対電流感度行列の各列に対して、最大値を有する非対角要素を識別して、この要素に値１を割り当て、更に他の全ての非対角要素に値０を割り当てて、各列に正確に２つの非ゼロ要素を有する正方「２値電流感度行列」を得るサブステップと、
ｃ．その成分が、ステップＩＶにおいて得られた前記行列の各列及び各行それぞれにおける要素の総和に等しい、第１のベクトル（Ｃ）及び第２のベクトル（Ｒ）を計算し、その成分が前記第２のベクトル（Ｒ）の各成分を前記第１のベクトル（Ｃ）の対応する成分で除算することによって得られる第３のベクトル（Ｓ）を更に計算することによって、ステップＩＶにおいて得られた前記行列のうち前記第１の計測ユニットに対応する行及び列の位置を決定するサブステップであって、ステップＩＶにおいて得られた前記行列における前記第１の計測ユニットに対応する前記行及び列の位置が、最大値を有する前記第３のベクトル（Ｓ）の成分を識別することにより決定できるようになるサブステップと、
ｄ．ステップＶで得られた前記２値電流感度行列から前記第１の計測ユニットに対応する前記列を削除して、各特定の分岐（スラック分岐を除く）が関連する２つのノードを示す２つの非ゼロ要素を各列に有する前記ネットワークの接続行列を得るサブステップと、
を含む、請求項１又は２に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項６】
前記ステップＩＩＩの最尤推定は、多重線形回帰の形態で実装される、請求項１から５のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項７】
前記ステップＩＩＩの最尤推定は、連続する時間ステップに対応する２つの誤差項間の相関が、０．７から０．３の間の区間に含まれ、好ましくは０．６から０．４の間の区間に含まれ、最も好ましくは約０．５に等しいこと、及び非連続時間ステップに対応する２つの誤差項間の前記相関が、０．３から０．３の間の区間に含まれ、好ましくは０．２から０．２の間の区間に含まれ、最も好ましくは約０．０であることを仮定して実行される、請求項１から６のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項８】
移動体通信事業者によって提供される既存の商用ネットワークが、前記通信ネットワークとして作用する、請求項１から７のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項９】
前記計測ユニットは、前記通信ネットワーク経由でのネットワークタイムプロトコル（ＮＴＰ）を用いて同期される、請求項１から８のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項１０】
前記計測ユニットは各々、コントローラ及びバッファを備え、前記ステップＩは、前記計測ユニットによって分散方式で一体化して実行される、請求項１から９のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項１１】
前記ステップＩが完了した後、前記処理ユニットは、前記通信ネットワークにアクセスして、前記計測ユニットから前記タイムスタンプ付き電流値をダウンロードする、請求項１０に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項１２】
前記計測ユニットは各々、コントローラ及びワーキングメモリを備え、前記計測ユニットのうちの１つは、前記処理ユニットとして作用する、請求項１から１１のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項１３】
前記電流測定値は、ＡＣ電力の少なくとも１つの周期にわたって測定された平均値である、請求項１から１２のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【請求項１４】
前記電力ネットワークは、アイランドモードで稼働し、前記電力ネットワークは、主グリッドから切り離される、請求項１から１３のいずれか１項に記載の電力ネットワーク（１）のトポロジ識別方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、概して、電力ネットワークの分野に関し、より具体的には、ネットワークパラメータ及びネットワークモデルの情報を用いずに電力ネットワークのトポロジを推定するための方法に関し、この電力ネットワークは、ノードのセット及び分岐のセットを含み、各分岐が１つのノードを別のノードに接続するように構成され、電力ネットワークは更に、計測ユニットと、通信ネットワークに接続された処理ユニットとを含む監視インフラストラクチャを備えており、計測ユニットは更に、この通信ネットワークに接続されて、処理ユニットとの間のデータ送信が可能になる。
【背景技術】
【０００２】
ネットワークのトポロジについての情報は、全ての電力システムの調査並びに技術的に安全で経済的に最適なグリッド運用にとって非常に重要な情報である。実際に、物理的グリッド構造についての情報は、ネットワーク監視及び状態推定、並びにグリッドの最適な制御及び管理などの、全ての電力システム調査に基本的な要素である。グリッドトポロジの識別は、信頼性の観点から及び多くのスマートグリッドアプリケーションにとって非常に重要である。一般に、配電システムのオペレータは、特に、低電圧ネットワークの場合に、ネットワークトポロジについての正確で最新の情報を有していない場合がある。更に、配電グリッドは、通常、放射状に稼働するが、１つの放射状トポロジから別のトポロジへの切り替えは、架線作業員によって行われる場合があり、場合により、スイッチ／ブレーカの状態は、この情報が時間内にシステムオペレータに到達することなく変更される可能性がある。ネットワークトポロジが変更された場合に、システムオペレータは、トポロジ変更後、スイッチ／ブレーカの状態を検証する必要がある場合がある。従って、グリッドトポロジについての更新された情報は、技術的に安全で経済的に最適なグリッド運用のための重要な情報である。更に、電気量の測定に基づくグリッドトポロジの自動的識別は、多くのスマートグリッドアプリケーション、特にプラグアンドプレイ機能を有するこれらのアプリケーションにとって非常に重要である。
【０００３】
グリッドトポロジの識別情報についての論題は、広く調査されている。利用可能な文献では、配電ネットワークは、基本的に、ライン（又は分岐）で接続されたバス（又はノード）のセットとしてモデル化される。バスは、配電ネットワーク内の分岐点、変電所、キャビネット、発電機、又は負荷を表し、バス間に延びるライン（又は分岐）は、１つのバスから別のバスに電力を搬送する配線を表す。研究者は、グリッドトポロジを識別するための様々な手法を使用している。これらの手法は、
ｉ）電力線通信（ＰＬＣ）信号［１］［２］［３］［４］［５］
ｉｉ）フェーザ測定ユニット（ＰＭＵ）測定［６］［７］［８］［９］
ｉｉｉ）計測インフラストラクチャからの測定データ［１０］［１１］［１２］［１３］［１４］［１５］［１６］［１７］［１８］［１９］［２０］
ｉｖ）ネットワークモデルの情報を使用したスイッチ状態の推定［２１］［２２］［２３］［２４］［２９］
ｖ）最終顧客のエネルギーメーターデータ［２５］［２６］［２７］、及び
ｖｉ）エネルギー価格信号［２８］
を含む。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】米国特許出願公開第2010/0007219号明細書
【特許文献２】米国特許出願公開第2015/0241482号明細書
【特許文献３】米国特許出願公開第2016/0109491号明細書
【特許文献４】米国特許出願公開第2015/0153153号明細書
【特許文献５】米国特許出願公開第2015/0052088号明細書
【特許文献６】米国特許出願公開第2007/0005277号明細書
【特許文献７】欧州特許出願公開第3065250号明細書
【特許文献８】米国特許出願公開第2014/0039818号明細書
【特許文献９】国際公開第2013/026675号
【特許文献１０】中国特許第104600699号明細書
【特許文献１１】中国特許出願公開第201510973479号明細書
【特許文献１２】米国特許出願公開第2012/0089384号明細書
【特許文献１３】米国特許出願公開第2013/0035885号明細書
【特許文献１４】カナダ国特許出願公開第2689776号明細書
【特許文献１５】欧州特許出願公開第2458340号明細書
【特許文献１６】欧州特許出願公開第3086093号明細書
【特許文献１７】米国特許出願公開第2014/0032144号明細書
【非特許文献】
【０００５】
【非特許文献１】T. Erseghe and S. Tomasin, ”Plug and Play Topology Estimation via Powerline Communications for Smart Micro Grids,” IEEE Transactions on Signal Processing , vol. 61, no. 13, pp. 3368 - 3377, 2013.
【非特許文献２】L. Lampe and M. O. Ahmed, ”Power Grid Topology Inference Using Power Line Communications,” in IEEE International Conference on Smart Grid Communications, Vancouver, 2013.
【非特許文献３】M. O. Ahmed and L. Lampe, ”Power Line Communications for Low-Voltage Power Grid Tomography,” IEEE Transactions on Communications , vol. 61, no. 12, pp. 5163 - 5175, 2013.
【非特許文献４】J. Bamberger, C. Hoffmann, M. Metzger and C. Otte, ”Apparatus, method, and computer software for detection of topology changes in electrical networks”. WO Patent WO 2012031613 A1, 15 March 2012.
【非特許文献５】S. Bolognani, N. Bof, D. Michelotti, R. Muraro and L. Schenato, ”Identification of power distribution network topology via voltage correlation analysis,” in 52nd IEEE Conference on Decision and Control, Florence, 2013.
【非特許文献６】B. Nicoletta, M. Davide and M. Riccardo, ”Topology Identification of Smart Microgrids,” 2013.
【非特許文献７】M. H. Wen, R. Arghandeh, A. von Meier, K. Poolla and V. O. Li, ”Phase Identification in Distribution Networks with micro-synchrophasors,” in IEEE Power & Energy Society General Meeting, Denver, 2015.
【非特許文献８】D. Deka, S. Backhaus and M. Chertkov, ”Estimating Distribution Grid Topologies: A Graphical Learning based Approach,” in Power Systems Computation Conference, Genova, 2016.
【非特許文献９】X. Li, H. V. Poor and A. Scaglione, ”Blind Topology Identification for Power Systems,” in IEEE International Conference on Smart Grid Communications, Vancouver, 2013.
【非特許文献１０】V. Arya, T. S. Jayram, S. Pal and S. Kalyanaraman, ”Inferring connectivity model from meter measurements in distribution networks,” E-Energy, pp. 173-182, 2013.
【非特許文献１１】S. Jayadev, N. Bhatt, R. Pasumarthy and A. Rajeswaran, ”Identifying Topology of Power Distribution Networks Based on Smart Meter Data,” in https://arxiv.org/pdf/1609.02678, 2016.
【非特許文献１２】V. Kekatos, G. B. Giannakis and R. Baldick, ”Online Energy Price Matrix Factorization for Power Grid Topology Tracking,” IEEE Transactions on Smart Grid , vol. 7, no. 3, pp. 1239 - 1248, 2016.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上記に列記した様々な手法の中で、計測インフラストラクチャからの測定データを使用するものが、本発明に最も適している。計測データからのグリッドトポロジの識別は、時系列の電圧変化の間の相関分析［１０］［１１］、電圧ベースの顧客クラスタリング［１６］、電圧変動に対する電力変動の相関分析［１５］、同期されたエネルギー測定値間の相関係数［１７］、電流振幅の同期された測定値と電流の位相角の同期された測定値との間の相関分析［１８］、最適化ベースの手法［１２］［１３］［１９］、及びエネルギー測定値の主成分分析［１４］を使用した文献で研究されている。これらの方法の主な欠点は、
−ノイズの非ランダム性が考慮されていないこと、
−結果が、計測ユニットを接続するラインの長さによって不正確になる場合があること ［１０］［１１］［１５］、
−解決するために計算集約的であること［１２］［１３］［１９］
−ネットワーク内の全ての負荷を計測する必要があること［１４］［１９］、及び
−エネルギー損失が、推定におけるノイズ及び誤差の原因であること［１４］［１９］
であり、欠点は更に、以下のとおりである。
−［１７］において、監視デバイスのセットの同期されたエネルギー測定値間の相関係数が、監視デバイスの各ペア間の相互関係を識別するのに使用される。規定された閾値よりも大きい相関係数を有する監視デバイスが、反復手順を通じて、互いに直接リンクされているとみなされる。
−［１８］における方法は、電流振幅及び電流の位相角（すなわち、複素数値）の同期された測定に基づいている。この方法では、ノードに接続された全ての分岐の電流が測定されることが望ましい。電流測定値間の相関が、測定値と最適にマッチするネットワーク構成を得るのに使用される。この方法は、回帰分析及び最適化問題の解決に基づく反復処理を含む。この方法の主な欠点は、ｉ）全ての分岐電流が必要であり、このことが、測定デバイスを全てのノード及び分岐に使用する必要があることを意味すること、ｉｉ）欠落した分岐の電流を決定するのにキルヒホッフの電流法則が使用され、このことが、キルヒホッフの電流法則が使用される場合に常に、ノード分岐の接続性についての情報が既知であることを前提としていることを意味すること、及びｉｉ）調整すべき幾つかの回帰分析パラメータが存在することである。
【課題を解決するための手段】
【０００７】
従って、本発明の目的は、従来技術の上述の欠点及び問題を軽減することである。本発明は、添付の請求項１に従って配電ネットワークのトポロジを識別するための方法を提供することによってこの目的を達成する。
【０００８】
本明細書では、配電ネットワークのトポロジは、このネットワークの物理ノード間の接続性として規定される。ネットワークの接続性又はトポロジは、グラフ形式で表すことができ、グラフの各ノードは、物理ノードのうちの１つを表し、グラフの各エッジは、電力線（又は分岐）を表す。接続行列と呼ばれる特定のタイプの行列が、一般に、電力ネットワークのグラフを表すのに使用される。接続行列は、トポロジの観点からネットワークを記述するため、ネットワークのトポロジを認識することは、接続行列を認識することを意味し、逆もまた同様である。グラフの接続行列（Ａ）は、ノード上のエッジの関連性を表し、無向グラフの場合、以下のように規定される。
【数１】
2020533926.tif 000003
無向グラフでは、各エッジは、異なるノードに関連する２つの端部を有する。正確に２つの「１」が接続行列（Ａ）の各列に存在することが理解されるであろう。言い換えると、接続行列の各列の総和は、２に等しい。
【０００９】
本発明により、電力ネットワークは、少なくとも１つの分岐を通ってそれぞれのノードに流入又はそれから流出する電流を測定するように異なる物理ノードに配置された計測ユニットを備える。更に、ネットワークの接続性の決定は、特定のノードに接続された分岐における測定された電流の変動に最も感度の高い計測ユニットが、該特定のノードのすぐ上流の物理ノードに配置された計測ユニットであるという一般的仮定に基づく。
【００１０】
ネットワークの接続性を決定する前に、本発明の方法は、少なくとも１つの分岐を通って特定のノードに流入又はそれから流出する電流の変動を、少なくとも１つの分岐を通って別のノードに流入又はそれから流出する電流の付随する（ｃｏｎｃｏｍｉｔａｎｔ）変動に関連付ける電流感度係数を統計的に推定するタスクを含む。これらの電流感度係数は、以下に与えられる偏導関数の値の推定値として解釈することができ、
【数２】
2020533926.tif 000004
ここで、ｉ，ｊ∈｛１，．．．，Ｎ｝は、ネットワークの異なる物理ノードに配置された計測ユニットを識別し、Ｉi及びＩjは、少なくとも１つの分岐を通って２つの異なるノードに流入又はそれから流出する電流の強度を表す。
【００１１】
電流感度係数は、複数の計測ユニットによって提供されるタイムスタンプ付き電流強度測定値を使用して計算される。次に、電力ネットワークの接続性を表すグラフの接続行列が、この感度係数から推測できる。実際には、計測ユニットＢが、特定の物理ノードに位置し、計測ユニットＡが、該特定のノードのすぐ上流の物理ノードに位置する場合に、Ｂに対するＡの電流感度行列は、１に近い傾向があり、それは、該特定のノードに流入又はそれから流出する全ての電流が、そのすぐ上流のノードにも流れるためである。従って、計測ユニットＡは、計測ユニットＢの電流強度の変動を観測する。その一方、Ａに対するＢの電流感度係数は、０に近く、それは、計測ユニットＢが、計測ユニットＡの下流に位置し、上流ノードに流入又はそれから流出する電流のほとんどが、該特定のノードを通って流れないためである。計測ユニットＢは、上位の電流強度の変動のほとんどを観測することができない。従って、大きな電流感度係数（例えば、１に近い）は、直接又は間接的に上位に位置するノードを示す。
【００１２】
本発明の方法の主な利点の一部を以下に列挙する。
−ＰＭＵ測定値又はＰＬＣ信号は必要でない。
−トポロジは、短期間に対応する限られた量の測定データ、例えば、２００秒に対応する１００ｍｓの時間ステップを有する２０００個のデータポイントを使用して識別される。
−ネットワークのモデル及び／又はネットワークパラメータは必要でない。
−タイムスタンプ付き電流測定値のみが必要である。
−本アルゴリズムは、測定におけるノイズの存在に対して堅牢である。
−方程式系は、線形であるため、計算は、コンピュータ計算上効率的である（最適化ベースの手法と対照的）。
−本アルゴリズムは、測定デバイスが限られた数の物理ノードのみに展開されているために、ネットワークノードの一部に流入又はそれから流出する電流が監視されない場合においても、測定デバイス間の接続性（すなわち、トポロジ）を決定することができる。言い換えると、本アルゴリズムは、ネットワーク内のあらゆるところに測定デバイスを配置することを必要としない。
【００１３】
本発明により、幾つかの計測ユニットが、異なるネットワークノードに配置される。計測ユニットは、好ましくは６０ミリ秒から１０秒の間、最も好ましくは１００ミリ秒の時間間隔でタイムスタンプ付きの電流測定値を提供することができる。本発明の方法は、タイムスタンプ付き電流測定値のみを必要とすることに留意されたい。各計測ユニットは、或るノードに配置され、計測ユニットは、このノードに接続された少なくとも１つの分岐又は複数の分岐の電流を測定する。従って、分岐が接続されているノードのうちの１つは既知であり、トポロジ識別の問題は、この分岐が接続されている他のノードを見つけることである。
【００１４】
本発明の目的は、測定データのみを使用することによって、計測デバイスを接続する物理ネットワーク、すなわちネットワークトポロジを識別することである。本発明の方法は、全てのネットワークノードに計測デバイスが展開されるとは限らない場合に、効果的に実施することができる。添付の図２及び３は、両方とも、４つのバスと３つの給電線とを備える同じ例示的なネットワークを示している。図２は、４つの計測ユニットがこのネットワークに割り当てられた（物理ノードごとに１つ）場合の例示的な事例を示している。図示のように、ノード１にある計測デバイスは、変圧器を通って流れる電流を測定し、ノード２、３、及び４にある計測デバイスは、それぞれ分岐Ａ、Ｂ、Ｃを通って流れる電流を測定する。図３は、３つの計測ユニットのみがネットワークに割り当てられた場合の別の事例を示している。前の事例と同様に、ノード１にある計測デバイスは、変圧器を通って流れる電流を測定し、ノード２及び４にある計測デバイスは、それぞれ分岐Ａ及びＣを通って流れる電流を測定する。しかしながら、計測ユニットは、残りのノードに配置されていない。図２及び３は更に、各々、図示の計測ユニット展開に対応する接続行列を含む。図２の場合、ネットワーク接続行列は、４×３行列であり、ここで、４はノード数を表し、３はノード間の分岐の数を表す。図３の場合、ネットワーク接続行列は、計測ユニット間の物理的接続性を反映する３×２行列である。
【００１５】
発明の方法は、電流感度係数に基づく。計測ユニットによって測定された電流強度の変動（ΔＩi）は、上位ノードにある計測ユニットによって及び／又はノードが位置するループ内の計測ユニットによって観測できるため、電流強度の変動は、ネットワークトポロジを識別するための基準として選択される。例えば、図２及び３に示されているネットワークの場合、「ノード４」における電流強度の変動は更に、「ノード３」及び「ノード１」で観測できるが、この電流強度の変動は、「ノード２」では確認できない。同様に、「ノード２」における電流強度の変動は、「ノード１」でのみ測定でき、この電流強度の変動は、「ノード３」及び「ノード４」では観測できない。或るノードの上流又は下流の分岐の電流強度の変動を観測することは、ネットワークパラメータ（すなわち、コンダクタのタイプ、サイズ、及び長さ）と無関係であり、その一方、或るノードにおける電圧変動は、ネットワークパラメータに応じて近接ノードを通して測定できることに留意する価値がある。従って、グリッドトポロジを識別するのに電流強度の変動を使用することは、電圧変動よりも効率的な手法である。更に、特定のノードに関連するネットワークの分岐を流れる電力（分岐電力）は、Ｐ＝ＵＩｃｏｓφで与えられ、ここで、Ｉは、この分岐を通って流れる電流の強度であり、Ｕは、この分岐が関連するノードの電圧であり、φは、該電圧と該電流との間の位相角であることに言及すべきである。当業者であれば、分岐電力変動は、この分岐電力変動がノード電圧変動の影響も受けるにもかかわらず、分岐電流強度の変動と同様に振る舞うことができることを理解するであろう。従って、電流強度の変動の利用は、分岐電力変動を使用するよりも好ましい。
【００１６】
本発明の方法はまた、ネットワークがアイランドモードで稼働している場合のトポロジを識別するのに使用することもできる。例えば、添付の図４は、「スイッチ」によって主グリッドから切り離され、アイランドモードで稼働するネットワークを示している。発電機「Ｇ１」が、アイランド稼働モードでネットワークを制御でき、この発電機が接続されたノード（「ノード２」）が、スラックノードとみなされると仮定する。この場合、「ノード１」にある計測ユニットは、全く電流を測定せず、スラックノード（「ノード２」）にある計測ユニットは、全ての計測ユニットの電流強度の変動を観測する。言い換えると、提案のアルゴリズムは、図４に示されているように、スラックノード及びネットワーク接続性行列を識別する。
【００１７】
分岐の電流強度の変動は、ノード電力変動に起因する。言い換えると、グリッドのあらゆるノードにおける電力変動は、このノードの上位にある分岐並びに同じループ内の分岐における電流強度の変動をもたらす。同様の方法が、分岐の電流強度の変動をノード電力変動に関連付ける電流電力感度係数に基づいてネットワークトポロジを識別するのに使用することができる。この代替電流電力方法は更に、本明細書の最後の例３で概説される。
【００１８】
それにもかかわらず、電流電力方法は、あらゆる計測ユニットが測定ノードにおける有効電力を計算することを必要とする。このことは、あらゆる測定ノードにおける計測ユニットが、例えばロゴスキーコイル（Ｒｏｇｏｗｓｋｉ  ｃｏｉｌ）又は変流器などの幾つかの電流測定デバイスを備えて、これらの測定ノード（すなわち、給電線）に接続されたあらゆる分岐の電流を測定することを意味する。図１０は、ネットワークと、ノード電力変動を計算するのに必要な電流測定値とを示している。例えば、分岐Ｂ及びＣは、両方ともノード３に関連する。従って、分岐Ｂを通ってノード３に流入又はそれから流出する電流の測定値
2020533926.tif 000005
及び分岐Ｃを通ってノード３に流入又はそれから流出する電流の測定値
2020533926.tif 000006
は、両方とも、「ノード３」における電力変動ΔＰ3を計算するのに必要である。更に、「ノード１」における電力変動ΔＰ1の計算は、電流の情報
2020533926.tif 000007
2020533926.tif 000008
及び
2020533926.tif 000009
を必要とする。この場合、多数の電流測定値をより効果的に使用してノード間の物理的接続性が識別できることが理解されるであろう。しかしながら、この場合、ネットワークトポロジは、電流電力感度係数を使用するのではなく、電流測定値を比較することによってより容易に識別することができる。このことは、分岐の両側の電流測定値が最も類似した特性を有するという事実に起因する。例えば、ノード１との接続点で分岐Ｂを通る電流の測定値（
2020533926.tif 000010
）は、ノード３との接続点で分岐Ｂを通る電流の測定値
2020533926.tif 000011
）と最も類似している。従って、以下の詳細な説明は、相互電流感度係数に基づく方法を開示することに当てられる。本方法は、限られた数の電流測定値を用いた実際の用途に最も有効である。
【００１９】
本発明の他の特徴及び利点は、単に非限定的な例として与えられ添付図面を参照して行われる以下の説明を読むことで明らかになるであろう。
【図面の簡単な説明】
【００２０】
【図１】本発明の方法の特定の実装形態を説明するのに使用される例示的な配電ネットワークの概略図である。
【図２】計測デバイスを同じ例示的な電気ネットワークに割り当てて、計測デバイスを結果として生じる異なる接続行列に割り当てる２つの代替例を示す。
【図３】計測デバイスを同じ例示的な電気ネットワークに割り当てて、計測デバイスを結果として生じる異なる接続行列に割り当てる２つの代替例を示す。
【図４】アイランドモードで動作する電気ネットワーク並びに対応する接続行列を示す。
【図５】電力ネットワークのトポロジを推定するための本発明の方法の第１の例示的な実装形態を示すフローチャートである。
【図６】本発明の方法の第２の特定の実装形態を示すフローチャートである。
【図７】本方法の第３の特定の実装形態を示すフローチャートである。
【図８】本発明の方法を用いてネットワークのトポロジを推定するのに使用される測定電流データの第１の例を示す。
【図９】本発明の方法を用いてネットワークのトポロジを推定するのに使用される測定電流データの第２の例を示す。
【図１０】図２及び図３に既に示されており電気ネットワークに関するノード電力変動を計算するのに使用できる計測デバイスの特定の割り当て例を示す。
【発明を実施するための形態】
【００２１】
本発明の主題は、電力ネットワークのトポロジを推定するための方法である。従って、本発明の適用分野は、電力ネットワークの分野であるので、例示的なネットワークが、最初に説明される。その後、本方法が動作できる実際の方法について説明する。
【００２２】
図１は、５７個の住居ブロックと、９つの農業施設とで構成され、合計８８人の顧客に供給を行う例示的な低電圧放射状配電ネットワーク（１で参照）の概略図である。低電圧ネットワーク１（２３０／４００ボルト、５０Ｈｚ）は、変電所変圧器によって中電圧ネットワーク３にリンクされる。図１では、変電所変圧器は、インピーダンスＺｃｃと組み合わされた理想変換器（５で参照）として表されており、このインピーダンスは、理想変換器５の出力とネットワーク１の残り部分との間に挿入される。以下の表１は、この特定の例での変電所変圧器に関する可能性のある特性の概念を与えることが意図されている。
【００２３】
表１
2020533926.tif 000012
【００２４】
変電所変圧器は、回路遮断器９及び第１のバスＮ１を介してネットワーク１に接続される。図示の例のネットワークでは、幾つかの給電線が、バスＮ１から分岐する。これらの給電線のうちの１つ（Ｌ１で参照）は、５つの住居ブロック及び１つの農業施設のサブセットを低電圧ネットワークにリンクするように構成される。残りの５２個の住居ブロック及び８つの農業施設は、図１に明示されていない（しかしながら、１１で参照される単一の矢印で全体として表されている）他の給電線によってバスＮ１にリンクできることを理解されたい。
【００２５】
給電線Ｌ１は、バスＮ１を第２のバス（Ｎ２で参照）に接続する。図１から理解できるように、３つの住居ブロック及び１つの農業施設が、バスＮ２に接続される。更に、給電線Ｌ２は、バスＮ２を第３のバス（Ｎ３で参照）に接続する。２つの住居ブロックがバスＮ３に接続される。表２（以下）は、この特定の例で使用される給電線Ｌ１及びＬ２に関する可能性のある特性の概念を与えることが意図されている。
【００２６】
表２
2020533926.tif 000013
【００２７】
引き続き図１を参照すると、ネットワーク１は更に、３つの分散型発電所を備えることが理解できる。第１の発電所（Ｇ１で参照）は、バスＮ２に接続された太陽光発電所であり、第２の発電所（Ｇ２で参照）は、バスＮ３に接続された太陽光発電所であり、第３の発電所は、バスＮ１にリンクされたディーゼル発電機である。更に詳細に説明するように、第３の発電所は、電力ネットワーク１がアイランドモードで動作している場合に電圧基準発電機として作用するように構成される。図１では、ディーゼル発電機は、インピーダンスＸｄと組み合わされた理想的な発電機（Ｇ３で参照）として表されており、このインピーダンスは、理想的な発電機Ｇ３の出力とネットワーク１の残りの部分との間に挿入される。表３Ａ及び表３Ｂ（以下）は、この特定の例で使用される３つの分散型発電所に関する可能性のある特性の概念を与えることが意図されている。
【００２８】
表３Ａ
2020533926.tif 000014
【００２９】
表３Ｂ
2020533926.tif 000015
【００３０】
この例によれば、太陽光発電所Ｇ１及びＧ２は、最大２２６ｋＶＡの電力を供給することが確認できる。図１はまた、ネットワーク１のバスＮ１に接続されたバッテリパック（１５で参照）も示している。３つの分散型発電所とバッテリーパック１５と回路遮断器９とを組み合わせた存在は、低電圧ネットワーク１を一時的にアイランド化する可能性を与える。以下の表４は、この特定の例で使用されるバッテリーパック１５に関する可能性のある特性の概念を与えることが意図されている。
【００３１】
表４
2020533926.tif 000016
【００３２】
本発明の方法が実施される物理的環境は更に、電力ネットワークに加えて監視インフラストラクチャを備える必要がある。本発明により、監視インフラストラクチャは、ネットワークのノードの選択において設けられた計測ユニットを備える。各計測ユニットは、特定のノードに位置し、このユニットは、給電線を通ってこの特定のノードに流入又はそれから流出する電流の強度を測定するように構成される。本発明によれば、ネットワーク内の同じノードに配置された幾つかの計測ユニットが存在でき、これらの計測ユニットの各々は、同じノードに関連する異なるライン（又は分岐）を通って流れる電流を測定するように構成される（以下の本文では、少なくとも１つの計測ユニットを備えているネットワークのノードは、「測定ノード」と呼ばれる）。本明細書で説明する例示的な実装形態により、測定ノードごとにただ１つの計測ユニットが存在する（すなわち、より正確には、ネットワークの位相ごとにただ１つの計測ユニットが存在する）。実際には、前述のように、図１に示されている例示的な低電圧電力ネットワーク１は、三相電力ネットワークである。電流を測定することに関する限り、三相電力ネットワークは、３つの単相電力ネットワークのセットとみなすことができる。このような場合、本発明の好ましい実装形態は、電流が、三相のそれぞれ１つに対して別々に測定されることを提供する。この測定は、３倍の数の計測ユニットを使用すること、又は代替的に、３つの異なる位相を別々に測定するように設計された計測ユニットを使用することによって行うことができる。
【００３３】
図１は、７つの異なる測定ノード（Ｍ１からＭ７で参照）の位置を示している。しかしながら、本発明によれば、任意の数の測定ノード、場合により、わずか２つの測定ノードが存在できることを理解されたい。更に、図１に示されている特定のネットワークに関して、詳細に示されていないネットワーク１の残りの部分は、場合により、追加の測定ノードを備えることができることを理解されたい。ノードＭ１からＭ７の計測ユニットは各々、少なくとも１つの電流、好ましくは幾つかの異なる電流を局所的に測定するように構成される。再度図１を参照すると、第１の測定ノードＭ１が、変電所変圧器をバスＮ１に接続していることが理解できる。第２の測定ノードＭ２は、バッテリーパック１５をバスＮ１に接続し、第３の測定ノードＭ３は、ＰＶシステムＧ２をバスＮ３に接続し、第４の測定ノードＭ４は、ＰＶシステムＧ１をバスＮ２に接続し、第５の測定ノードＭ５は、ディーゼル発電機をバスＮ１に接続し、第６の測定ノードＭ６は、給電線Ｌ２をバスＮ３に接続する。最後に、第７の測定ノードＭ７は、給電線Ｌ１をバスＮ２に接続する。
【００３４】
本発明により、監視インフラストラクチャは更に、通信ネットワークを備え、計測ユニットは、この通信ネットワークに接続されて、処理ユニット７との間のデータ送信が可能になる。図１の非常に概略的な図では、処理ユニット７は、ネットワーク１から或る距離に配置されたコンピュータの形態で表される。しかしながら、処理ユニットは、測定ノードのうちの１つに位置することができることが理解されるであろう。実際には、監視インフラストラクチャの好ましい実施形態によれば、処理ユニットは、計測ユニットのうちの１つの一部分を形成する。図示の例によれば、通信ネットワークは、専用の伝送ネットワークではなく、移動体通信事業者によって提供される商用ＧＳＭネットワークである。しかしながら、代替実装形態により、監視インフラストラクチャ用の通信ネットワークは、当業者が適切とみなす任意のタイプのものであり得ることが理解されるであろう。
【００３５】
図５は、電力ネットワークのトポロジを推定するための本発明の方法の第１の例示的な実装形態を示すフローチャートである。図５のとりわけ詳細でないフローチャートは、４つのボックスで構成される。最初のボックス（０１で参照）は、一般に、少なくとも１つの分岐を通って電力ネットワークの複数のノードのそれぞれ１つに流入又はそれから流出する電流の強度に関する一連の値を取得することにあるタスクを表す。この目的のために、本発明の方法は、ネットワーク内の幾つかの位置で或る期間τにわたって繰り返し分岐電流を測定するように構成された監視インフラストラクチャを使用する。本発明により、電力ネットワークは、ＡＣ電力ネットワークであり、ほとんどの実装形態によれば、電流の測定値は、瞬時値でなく、ＡＣ電力の少なくとも半周期にわたって、好ましくはＡＣ電力の２周期から１０周期にわたって、最も好ましくはＡＣ電力の３周期にわたって（すなわち、５０Ｈｚ  ＡＣ電力ネットワークの場合、６０ｍｓの間）測定された平均値（好ましくは信号の基本周波数に基づくｒｍｓ値）である。本発明の方法は、異なる計測ユニットの測定値が高度に同期されることを必要としない。しかしながら、本方法は、異なる測定ノードにおける計測ユニットが、ほぼ同時に得られた測定値を提供することを必要とし、すなわち、言い換えると、本方法は、異なる測定ノードにおける測定が、得られた値を付随するものとして後で処理できるようにするのに十分に近い時間に行われることを必要とする。
【００３６】
本明細書で説明する本発明の実装形態により、ネットワーク内の異なる計測ユニットは、監視インフラストラクチャ用の通信ネットワークとして作用するＧＳＭネットワーク経由でネットワークタイムプロトコル（ＮＴＰ）を用いて同期される。ＮＴＰの利点は、実装するのが容易であり、ほぼいたるところで容易に利用できることである。ＮＴＰの公知の欠点は、非常に正確というわけではないことである。しかしながら、予期できるものに反して、経験により、ＮＴＰによって提供される同期は、本発明の方法が満足できる結果をもたらすのに十分好適であることが示されている。しかしながら、ＮＴＰは、本発明の方法と共に使用できる唯一の同期方法ではないことを理解されたい。具体的には、非常に費用のかかる実装形態によれば、計測ユニットは、共通の時間基準又はＧＰＳ同期との常設リンクを有するＰＭＵとすることができる。
【００３７】
本発明により、異なる計測ユニットによって行われる電流の測定は、上述した程度に同期される。この例により、計測ユニットは、好ましくは所与の時間ウィンドウ内で一定間隔で、電流を繰り返し測定する。連続する測定回数は、好ましくは２００から５０００回の間の測定に含まれ、好ましくは１０００から３０００回の間の測定に含まれ、例えば、２０００回の測定である。しかしながら、最適な測定回数は、測定ノード数の関数として増加する傾向があることを理解されたい。一方で、最適な測定回数は、計測ユニットによって提供される測定値の精度が高まるとともに、並びに計測ユニット間の同期の精度が高まるとともに、減少する傾向がある。
【００３８】
この例などでは、計測ユニットによって測定される値は、瞬時値でなく、ＡＣ電力の少なくとも半周期にわたって測定される平均値であり、連続する測定間の最小時間間隔は、ＡＣ電力の幾つかの周期に等しい必要がある。実際には、第１の例示的な実装形態によれば、連続する測定を隔てる時間間隔の長さは、好ましくは６０ｍｓから３秒の間であり、最も好ましくは６０ｍｓから１秒の間である。
【００３９】
図５のフローチャートにおける２番目のボックス（０２で参照）は、最初に、少なくとも１つの分岐を通って、ノードの各々に流入又はそれから流出する電流を測定するように異なる物理ノードに配置された計測ユニットのそれぞれ１つによって取得されたデータに基づいて、測定された電流強度の変動を計算するタスクと、更に、計測ユニットのそれぞれ１つによって測定された電流の強度の変動の、他の全ての計測ユニットによって測定された電流の強度の付随する変動に関するテーブルを編集するタスクとを表している。電流の強度の変動は、電流の各測定値から同じ変数の以前の値を減算することによって、計算することができる。言い換えると、同じ計測ユニットによる２つの測定値が、時間ｔ及びｔ＋Δｔで利用可能である場合に、各測定ノードとの接続点での分岐電流の変動
2020533926.tif 000017
は、
2020533926.tif 000018
として計算され、ここで、ｉ∈｛１，．．．，Ｎ｝は、ｉ番目の測定ノードに配置された計測ユニットを指定する。更に、本明細書では、測定値に対応する量が、ティルデ（例えば、
2020533926.tif 000019
）で示されることに留意されたい。
【００４０】
前述のように、本発明の非常に費用のかかる実装形態によれば、計測ユニットは、共通の時間基準又はＧＰＳ同期との常設リンクを有するＰＭＵとすることができる。この場合、電流の振幅及び位相の両方が測定される。電流の位相に関する情報が更に利用可能である場合には、必要な連続測定の回数は、電流の絶対値及び位相の両方を考慮することにより減少することが可能であり得る。実際には、この場合、測定された電流Ｉi（ｔ）は、複素数として処理でき、電流の２つの連続した測定値の間の差分が、同様に複素数として処理することができる。この場合、電流の変動
2020533926.tif 000020
は、好ましくは電流の２つの連続する測定値の差分に対応する複素数の絶対値として、すなわち、言い換えると、２つの連続する電流フェーザ間の差分の大きさとして計算される。
【００４１】
ここで本発明の第１の例示的実装形態に戻ると、処理ユニットは、電流強度の変動を計算するために、最初に通信ネットワークにアクセスして、異なる計測ユニットのバッファから電流
2020533926.tif 000021
に関するタイムスタンプ付き値をダウンロードする。次に、処理ユニットは、ダウンロードされた各電流値から、直前のタイムスタンプを保持する同じ変数の値を減算することによって、測定された電流強度の変動を計算する。特に、時間ｔ∈｛ｔ1，．．．，ｔm｝は、異なる計測ユニットによって提供されるタイムスタンプを指すことに留意されたい。例えば、Ｉ1（ｔ1）及びＩN（ｔ1）は、異なる計測ユニットからの測定値から計算され、第１の例示的な実装形態により、これらの計測ユニットのそれぞれのクロックは、ＮＴＰを使用して同期されるため、時間ｔでの測定値は、時間ｔ±標準ＮＴＰ同期誤差における測定値を意味するものとして理解されたい。
【００４２】
次に、処理ユニットは、少なくとも１つの分岐を通って１つの特定の測定ノードに流入又はそれから流出する電流の測定された強度の各変動
2020533926.tif 000022
（ここで、ｉ∈｛１，．．．，Ｎ｝は、ｉ番目の測定ノードに配置された計測ユニットを指定する）を、他の全ての測定ノードに配置された計測ユニットによって同じ時間に（ここで、ｔ∈｛ｔ1，．．．，ｔm｝は、特定の測定時間又はタイムスタンプを表す）測定された電流の強度の変動と関連付ける。表５（下記）で例示されるように、結果は、表として表すことができ、この表の列は、異なる測定ノードｉに配置された異なる計測ユニットに対応し、この表の行は、連続する測定時間に対応する。これらの測定時間は、所与の時間ウィンドウτ＝［ｔ1，ｔm］をカバーする。本発明により、ｍ＞２Ｎであり、好ましくはｍ＞＞Ｎである。
【００４３】
表５
2020533926.tif 000023
【００４４】
図５のフローチャートにおける３番目のボックス（０３で参照）は、特定の計測ユニットによって測定された電流強度の変動を、他の計測ユニットによって測定された電流強度の変動に関連付ける電流感度係数を統計的に推定するタスクを表す。表５のデータから得られる電流感度係数のセットは、好ましくは最尤推定（ＭＬＥ）法を用いて得られる。これらの電流感度係数はグループ化されて、電流感度係数行列が形成できる。１つの計測ユニットによって測定された電流強度の変動を、他の計測ユニットによって測定された付随する電流強度の変動に関連付ける電流感度係数は、以下に与えられる偏導関数の値の推定値として解釈することができる。
【数３】
2020533926.tif 000024
【００４５】
最尤推定によると、他のユニットに対する各計測ユニットの電流感度係数は、以下の最適化問題又はその凸再公式化の結果として得ることができ、
【数４】
2020533926.tif 000025
（以下の条件）
【数５】
2020533926.tif 000026
ここで、ΔＩ（ｔ）iは、測定された電流強度の変動であり、
2020533926.tif 000027
は、ノイズのない推定電流強度の変動であり、
2020533926.tif 000028
である。
【００４６】
計測ユニットの、それ自体に関する電流感度係数は、１に等しい、すなわち、ＫＩＩi,j＝１である。分岐の電流感度係数が、その分岐自体を含む全ての分岐に対して計算される場合には、その分岐自体に関する係数は、１に等しくなり、他の分岐に関する係数は、無視できるようになり、それは、時系列時間データとそのデータ自体との間の相関が１に等しいためである。言い換えると、あらゆる計測ユニットは、その計測ユニット自体の測定値に関して他の計測ユニットよりも最高の感度係数を有する。或る計測ユニットの、他の計測ユニットに対する電流感度係数を決定するために、係数が計算される分岐の測定電流強度の変動は、推定においてΔＩ（ｔ）iから除外される必要がある。例えば、計測ユニットｊに関するＫＩＩi,jは、ΔＩ（ｔ）i（すなわち、ΔＩi≠j（ｔ））からｊ番目の列を除外することによって計算され、ＫＩＩi,jは、１に等しいとみなされる。
【００４７】
本方法の統計的な性質に起因して、個々の測定値は、その予測値からある程度逸脱する傾向がある。従って、各測定される電流強度の変動は、対応する推定電流強度の変動プラス／マイナス誤差項に等しい。すなわち、
【数６】
2020533926.tif 000029
であり、ここで、ωi（ｔ）は、誤差項である。
【００４８】
本発明によれば、最尤推定（ＭＬＥ）は、負の一次自己相関を考慮する。このことは、ＭＬＥが、誤差ωi（ｔ）と誤差ωi（ｔ＋Δｔ）との間に実質的な負の相関が存在すると仮定していることを意味し、ここで、ｔ及びｔ＋Δｔは、２つの連続する時間ステップである。本明細書において、「実質的な相関」という表現は、その大きさが、少なくとも０．３であり、好ましくは少なくとも０．４であり、最も好ましい場合に、ほぼ０．５に等しい相関を意味することを意図している。
【００４９】
発明の好ましい実装形態により、ＭＬＥは更に、２つの非連続時間ステップからの誤差間に実質的な相関が存在しないことを仮定する。「実質的な相関がない」という表現は、その大きさが、０．３より小さく、好ましくは０．２より小さく、最も好ましい場合に、ほぼ０．０に等しい相関を意味することを意図している。従って、２つの非連続時間ステップにおける誤差間の相関は、−０．３から０．３の間の区間に含まれ、好ましくは−０．２から０．２の間の区間に含まれ、この相関は、最も好ましい場合に、ほぼ０．０に等しい。連続する測定の回数は、ｍであるため、計測ユニットごとにｍ−１個の誤差項ωi（ｔ）が存在し、従って、（ｍ−１）×（ｍ−１）個の誤差相関項が存在する。
【００５０】
図５のフローチャートにおける４番目のボックス（０４で参照）は、前のステップで得られた電流感度係数行列からネットワーク接続性を導出するタスクを表す。接続性の決定は、少なくとも１つの分岐を通って特定のノードに流入又はそれから流出し、別の計測ユニットによって測定された電流の変動に最も感度の高い計測ユニットが、この特定のノードのすぐ上流のノードに配置されたものであるという一般的原理に基づく。言い換えると、配電ネットワークでは、大きな電流感度係数（例えば、１に近い）は、直接又は間接的に上流に位置するノードを示す。
【００５１】
図６は、図５のフローチャートで示される実装形態の特定の変形形態を示すフローチャートである。図示の変形形態によれば、多重尤度推定（Ｍｕｌｔｉｐｌｅ  Ｌｉｋｅｌｉｈｏｏｄ  Ｅｓｔｉｍａｔｉｏｎ）は、多重線形回帰（ｍｕｌｔｉｐｌｅ  ｌｉｎｅａｒ  ｒｅｇｒｅｓｓｉｏｎ）の形態で実装される。更に具体的に言うと、実装される特定のタイプの多重線形回帰は、「一般化最小二乗法」である。一般化最小二乗法は、測定ノードｉ∈｛１，．．．，Ｎ｝に配置された各計測ユニットに対して、以下の方程式を解くことを通じて分析的に電流感度係数を得ることを可能にし、
【数７】
2020533926.tif 000030
【数８】
2020533926.tif 000031
ここで、ΔＩ（ｔm）iは、時刻ｍでのｉ番目の分岐の電流強度の変動であり、Σm,mは、１次自己相関を用いて測定ノイズの影響を考慮するための相関行列である。
【００５２】
一般化最小二乗多重線形回帰法の結果は、誤差、すなわち、ωi（ｔ）が、既知の共分散行列を有する多変量正規分布に従う場合に、最尤推定（ＭＬＥ）と同じである。誤差相関行列Σm,mは、好ましくは処理ユニットに事前にロードされず、測定電流の変動についてのテーブル（表５）が生成された後にのみ、生成される（ボックス０２）。実際には、（ｍ−１）×（ｍ−１）誤差相関行列のサイズは、測定された電流の変動のテーブルの長さｍ−１によって定まる。従って、図６の変形形態は、図５に存在しない追加のボックス０２Ａを含む。ボックス０２Ａは、各測定ノードに配置された計測ユニットに関する誤差相関行列を生成するタスクを含む。ボックス０２の後のこの追加ボックスの存在は、１つの特定のタイムスタンプに関連するデータのセットが欠落している場合に、本方法を適合させることを可能にするという利点を有する。
【００５３】
この例では、任意の相関行列の場合と同様に、Ｎ（ｍ−１）×（ｍ−１）相関行列のそれぞれ１つの主対角要素のエントリは、全て、１に等しいように選択される。本発明によれば、下の第１の対角要素及びその上の第１の対角要素の両方のエントリは、全て、−０．７から−０．３の間に含まれ、最後に、他の全てのエントリは、−０．３から０．３の間に含まれる。この特定の例では、２つの非連続時間ステップ間の誤差の相関係数は、ゼロに等しく、２つの連続する時間ステップ間の誤差の相関係数は、−０．５であると仮定される。この場合、誤差相関行列は、以下に示されている三重対角行列に対応する。
【数９】
2020533926.tif 000032
【００５４】
図７は、図５のフローチャートで示されている実装形態の別の特定の変形形態を示すフローチャートである。図７の変形形態によれば、ＫＩＩi,jで与えられる電流感度係数行列からネットワーク接続性を推論するステップ（図５におけるステップ［０４］）は、２つの異なるステップ（［０４Ａ］及び［０４Ｂ］）に分解される。ステップ［０４Ａ］は、「スラックノードを見つけて、スラック電流に対応する列を絶対感度行列から削除するステップ」に当てられる。この部分に対してステップ［０４Ｂ］は、「ネットワーク接続行列を計算するステップ」に当てられる。図７に示されている変形形態は、計測ユニットがスラックノードに取り付けられていると仮定する。スラックノードは、本発明の方法が適用されるネットワークをより高い電圧ネットワークに接続する変圧器の低電圧側とすることができる。スラックノードとしての変圧器の場合、変圧器電流は、スラック分岐電流とみなされる。従って、スラック分岐は、スラックノードの上流にある。放射状ネットワークにおけるこの仮定は、ノード数（Ｎ）が、ネットワーク内の分岐の数に１を加えたものに等しいことを意味し、ここで、１は、スラック電流を測定する計測ユニットを表す。放射状ネットワークの場合、ネットワーク接続行列は、Ｎ×（Ｎ−１）行列であり、ここで、行の数は、スラックノードを含むノードの数を示し、列の数は、スラック分岐を除く行の数を示すことに留意されたい。
【００５５】
図７におけるボックス［０４Ａ］は、以下の手順で実装することができる。
ａ）電流感度行列の各係数ＫＩＩi,jをその絶対値で置き換えることにより、「絶対電流感度行列」を生成する。
ｂ）絶対電流感度行列の各列の要素を合計して、その結果を行ベクトルＣ1×Nに保存する。同様に、絶対電流感度行列の各行の要素を合計して、その結果をＲN×1の列ベクトルに保存する。次に、「スラックノードインジケータベクトル」ＳN×1を以下のように計算する。
【数１０】
2020533926.tif 000033
ｃ）最大値を有するスラックノードインジケータベクトルの成分を見つける。ｉで示されるこの成分のインデックスは、スラックノードにある計測ユニットを示す。絶対電流感度行列における対応する列は、スラック分岐を示す。この手順は、スラックノードの計測ユニットが、ネットワークの全ての分岐における電流強度の変動を観測するが、スラックノードにおける電流強度の変動が、ネットワーク全体を通して他の計測ユニットによって観測されないという事実に起因する。
【００５６】
図７におけるボックス［０４Ｂ］は、以下の手順で実装することができる。
ａ）絶対電流感度行列を取得し、この行列の、ｊ≠ｉである各列「ｉ」に対して、最大値を有する要素を識別して、この要素を１に等しく設定し、更に同じ列の要素
2020533926.tif 000034
の値を値１で置き換える。最後に、各列における他の要素の値を０に等しく設定する。これらの変換の結果は、以下の式、すなわち、
【数１１】
2020533926.tif 000035
に対応する行列
2020533926.tif 000036
を与える
（
2020533926.tif 000037
の各列は、２つの非ゼロ要素を有する）。
ｂ）Ｎ×Ｎ行例
2020533926.tif 000038
から、ステップ［０４Ａ］におけるスラック分岐のインデックスに対応するインデックス「ｊ」を有する列を削除する。結果として生じるＮ×（Ｎ−１）行例は、ネットワーク接続行列である。上述のように、各列は、２つの非ゼロ要素を有する。これらの２つの要素は、特定の列に対応する分岐が関連する２つのノードを示す。
【００５７】
上述のように、ステップ０４Ａ及び０４Ｂは、必ずしも連続して実行されるとは限らないことに留意することが重要である。更に、Ｎ×Ｎ行列から１つの列を削除してＮ×（Ｎ−１）を得るステップは、代替的に、電流感度係数を１及び０で置き換えて行列Ａを生成する前に、行うことができる。
【００５８】
例１
添付の図８は、４つの計測ユニットを有するネットワークの単線結線図（左側）と、すべての計測ユニットからのタイムスタンプ付き電流測定値（右側）とを示している。対応する電流感度係数行列は、上記で説明した一般化最小二乗多重線形回帰を用いて計算される。計算された電流感度係数行列は、次のとおりである。
【数１２】
2020533926.tif 000039
太字の数字で示されている対角要素は、各計測ユニットのそれ自体に関する係数を示しており、１に等しいことに留意されたい。絶対電流感度行列の「スラックノードインジケータベクトル」は、以下のとおりである。
【数１３】
2020533926.tif 000040
スラックノードインジケータベクトルは、「計測ユニット１」がスラックノード（すなわち、ノード１）に接続されていることを示す。従って、最初の列は、電流感度行列ＫＩＩから削除される。縮小された４×３行列が以下に与えられ、ここで、各列における最大感度係数は、計測ユニットのそれ自体に関する係数に関係なく、太字の数字で示されている。例えば、列１は、最大感度係数が、「計測ユニット２」と「計測ユニット１」との間にあり、これら２つのノードの間の物理的接続性が示される。同様の観察は、「計測ユニット３」と「計測ユニット２」との間の接続性に関して列２を考察すること、及び「計測ユニット４」と「計測ユニット１」との間の接続性に関して列３を考察することによって行うことができる。
【数１４】
2020533926.tif 000041
太字の数字で示されている要素は、ネットワーク接続行列（Ａi,j）における非ゼロ要素に対応する。
【数１５】
2020533926.tif 000042
【００５９】
例２
添付の図９は、図８と同じネットワークを示している。しかしながら、このネットワークは、３つの計測デバイスのみを備えている。この場合、計測ユニットに対して異なる番号付け順序が選択され、言い換えると、「計測ユニット１」、「計測ユニット２」、及び「計測ユニット３」が、それぞれ、ノード３、ノード１、及びノード４に配置されることに留意されたい。この場合、結果として生じる電流感度係数は、以下の行列で与えられる。この場合、結果として生じる電流感度係数は、以下の行列で与えられる。
【数１６】
2020533926.tif 000043
絶対電流感度行列の「スラックノードインジケータベクトル」は、
【数１７】
2020533926.tif 000044
である。この行列は、「計測ユニット２」がスラックノード（すなわち、ノード１）に接続されていることを示す。従って、２番目の列は、電流感度行列から削除される。縮小された３×２行列が以下に与えられ、ここで、各列における最大感度係数は、計測ユニットのそれ自体に関する係数に関係なく、太字の数字で示されている。
【数１８】
2020533926.tif 000045
従って、ネットワーク接続行列は、以下の上部で与えられる。
【数１９】
2020533926.tif 000046
列１は、「計測ユニット１」と「計測ユニット２」との間の物理的接続性を示し、列２は、「計測ユニット３」と「計測ユニット２」との間の物理的接続性を示している。興味深いことに、例２の結果は、全てのネットワークノードに計測デバイスが配置されているとは限らない場合、又は計測デバイスに対して番号付け順序が異なる場合でさえ、開発されたアルゴリズムが、複数の測定ノード間の物理的接続性を識別できることを示している。
【００６０】
例３
また、前述のように、ネットワークトポロジは、電流電力感度係数から推測することもできる。この代替電流電力方法について以下に概説する。本方法は、電力ネットワークが、スラックノード（Ｎ１１）及び他の計測ノード（Ｎ１２、．．．、Ｎ１４）を含む計測ノードのセット（Ｎ１１、．．．、Ｎ１４、図１０を参照）と、このスラックノードに接続されたスラック分岐（２０）と、これらの計測ノードを互いに接続するように構成された他の分岐のセット（Ａ、Ｂ、Ｃ）とを備え、電力ネットワークが更に、これらの測定ノードのそれぞれ１つの電圧を測定するための手段を含む監視インフラストラクチャを備えていることを必要とし、測定するための手段は、計測ユニット（Ｍ１１、．．．、Ｍ１７）を備え、これらの計測ユニットのうちの１つ（Ｍ１１）は、スラック分岐（２０）を通ってスラックノード（Ｎ１１）に流入又はそれから流出する電流を測定するように構成され、測定ノードのそれぞれ１つは更に、１つの測定ノードに関連する他の分岐（Ａ、Ｂ、Ｃ）のそれぞれ１つを通って測定ノードのうちの１つに流入又はそれから流出する電流を測定するための計測ユニット（Ｍ１２、．．．、Ｍ１７）を備えており、監視インフラストラクチャは、通信ネットワークに接続された処理ユニット（７）を備え、計測ユニットは、この処理ユニットとの間のデータ送信が可能になるようにこの通信ネットワークに接続され、電流電力方法は、
Ｉ．測定ノード（Ｎ１１、．．．、Ｎ１４）のそれぞれ１つの電圧を監視して、計測ユニット（Ｍ１１、．．．、Ｍ１７）に、時間ウィンドウ（τ）にわたって繰り返して電流強度値（Ｉ（ｔ））及び電流強度値のタイムスタンプｔ∈｛ｔ1，．．．，ｔm｝を同時に測定させるステップと、
ＩＩ．各測定ノード（Ｎ１１、．．．、Ｎ１４）に対してノード電力Ｐiを計算し、更に、各ノード電力値から以前のノード電力値を減算することによって、測定ノードのそれぞれ１つのノード電力の変動
2020533926.tif 000047
を計算するステップと、
ＩＩＩ．各電流強度値から、同じ計測ユニットによって測定された以前の電流強度値を減算することによって、測定ノードのそれぞれ１つに提供されて計測ユニットのうちの少なくとも１つによって測定された電流強度の変動ΔＩi（ｔ）を計算し、測定ノード（Ｎ１１、．．．、Ｎ１４）のそれぞれ１つにおける電流強度の変動ΔＩi（ｔ）の、全ての測定ノードにおけるノード電力の付随する変動
2020533926.tif 000048
に関する時系列順序付きテーブルを編集するステップと、
ＩＶ．実際の電流強度の変動（ΔＩi（ｔ））と最尤推定によって予測される変動との間の差異に対応する誤差項間の負の１次系列相関を考慮しながら、ステップＩＩＩの間に編集されたように、各測定ノード（Ｎ１１、．．．、Ｎ１４）において測定された電流強度の変動ΔＩi（ｔ）を、全ての測定ノードにおけるノード電力変動ΔＰj（ｔ）に関連付ける電流電力感度係数ＫＩＰi,jの最尤推定を実行して、電流電力感度係数行列ＫＩＰi,jを得るステップと、
Ｖ．ステップＩＶで推定された電流電力感度係数からネットワーク接続行列ＡN×(N-1)を計算するステップと、を含む。
【００６１】
電流電力方法の好ましい実装形態により、ステップＶは、以下のサブステップ、すなわち、
ａ．電流電力感度係数行列の各要素に、対応する要素の推定値の絶対値を割り当てることによって「絶対電流電力感度行列」を得るサブステップと、
ｂ．絶対電流電力感度行列の各列に対して、最大値を有する非対角要素を識別して、この要素に値１を割り当て、更に他の全ての非対角要素に値０を割り当てて、更に各対角要素に値１を割り当てることにより、各列に正確に２つの非ゼロ要素を有する正方「２値電流電力感度行列」を得るサブステップと、
ｃ．その成分が、ステップＩＶにおいて得られた行列の各列及び各行それぞれにおける要素の総和に等しい、第１のベクトル（Ｃ）及び第２のベクトル（Ｒ）を計算し、その成分が第２のベクトル（Ｒ）の各成分を第１のベクトル（Ｃ）の対応する成分で除算することによって得られる第３のベクトル（Ｓ）を更に計算することによって、ステップＩＶにおいて得られた行列のうちの第１の計測ユニットに対応する行及び列の位置を決定するサブステップであって、ステップＩＶにおいて得られた行列における第１の計測ユニットに対応する行及び列の位置が、最大値を有する第３のベクトル（Ｓ）の成分を識別することにより決定できるようになるサブステップと、
ｄ．ステップＶで得られた２値電流電力感度行列から第１の計測ユニットに対応する列を削除して、各特定の分岐（スラック分岐を除く）が関連する２つのノードを示す２つの非ゼロ要素を各列に有するネットワーク接続行列を得るサブステップと、を含む。
【００６２】
図１０は、図８及び９にも示されている同じネットワークの単線結線図を示している。
しかしながら、図８に示されているネットワークの各ノードは、単一の計測ユニットを備えているため、４つの計測ユニットがネットワーク全体に存在し、その一方、図１０に示されているネットワークには７つの計測ユニットが存在する。図１０における追加の計測ユニットが、ノード電力を計算するのに必要である。３つの計測ユニットがノード１に位置し、２つの計測ユニットがノード３に位置し、１つの計測ユニットがノード２及びノード４のそれぞれ１つに位置する。「ノードｉ」におけるノード電力（Ｐｉ）は、以下のように計算でき、
【数２０】
2020533926.tif 000049
ここで、Ｕｉはノードｉの電圧であり、Ｉkは、ノードｉに関連する分岐ｋを通って流れる電流であり、φkは、分岐ｋにおける電流と電圧との間の位相角である。
【００６３】
次のステップは、各ノードにおけるノード電力の変動
2020533926.tif 000050
並びに、少なくともネットワークの分岐の選択を通じた測定電流の付随する変動
2020533926.tif 000051
を計算して、選択された分岐のそれぞれ１つを通る電流の変動の、全てのノードにおける電力の付随する変動に関するテーブルを編集するためのものである。前述のように、この例では、７つの計測ユニットが、ネットワークを監視するのに使用される。
【００６４】
対応する電流電力感度行列の係数は、以下に与えられる偏導関数の値の推定値として解釈することができる。
【数２１】
2020533926.tif 000052
上記で説明したように、この係数は、上記で説明した一般化最小二乗多重線形回帰を用いて計算することができる。回帰分析は、７つの計測ユニット全てによって測定された電流の強度の変動値を予測することを可能にする。しかしながら、前の例では、任意の所与のノードに配置された計測ユニットは、単一より多くなかったため、パラメトリック重回帰分析は、選択された計測ユニットの各々が異なるノードに位置する４つの計測ユニットを選択することによって測定される電流の強度に影響を与える変動のみが予測されるように適用されるように選択される。本例では、各ノードに対して、ノードの上流に配置された分岐を通ってノードに流入又はそれから流出する電流を測定するように構成された計測ユニットが選択される。従って、計算された電流電力感度係数行列は、以下のとおりである。
【数２２】
2020533926.tif 000053
ｉ番目の行且つｊ番目の列における要素は、ｊ番目の列の電力変動に対するｉ番目の計測ユニットの電流強度の変動を示す。例えば、ＫＩＰ1,4＝３．７７０８及びＫＩＰ4,4＝３．９７０９の大きな値は、ノード４における電力変動が、計測ユニット１及び４によって測定される電流の強度に影響を与えることを示す。その一方、ノード４における電力変動は、観測されず、計測ユニット２及び３によって測定される電流の強度に実質的に影響を与えない（ＫＩＰ2,4＝０．０００３及びＫＩＰ3,4＝０．０００４）。
【００６５】
本発明の方法について例示的な実装形態によってより詳細に図示し説明したが、本発明は、開示した例に限定されるものではなく、当業者は、添付の特許請求の範囲によって定まる本発明の範囲から逸脱することなくこれらの例から様々な変更例及び／又は改良例を得ることができる。
【００６６】
１−参考文献
[1]  T. Erseghe and S. Tomasin, ”Plug and Play Topology Estimation via Powerline Communications for Smart Micro Grids,” IEEE Transactions on Signal Processing , vol. 61, no. 13, pp. 3368 - 3377, 2013. 
[2]  L. Lampe and M. O. Ahmed, ”Power Grid Topology Inference Using Power Line Communications,” in IEEE International Conference on Smart Grid Communications, Vancouver, 2013. 
[3]  M. O. Ahmed and L. Lampe, ”Power Line Communications for Low-Voltage Power Grid Tomography,” IEEE Transactions on Communications , vol. 61, no. 12, pp. 5163 - 5175, 2013. 
[4]  J. Bamberger, C. Hoffmann, M. Metzger and C. Otte, ”Apparatus, method, and computer software for detection of topology changes in electrical networks”. WO Patent WO 2012031613 A1, 15 March 2012.
[5]  E. G. de Buda and J. C. Kuurstra, ”System for Automatically Detecting Power System Configuration”. US Patent US 20100007219 A1, 14 January 2010.
[6]  S. Bolognani, N. Bof, D. Michelotti, R. Muraro and L. Schenato, ”Identification of power distribution network topology via voltage correlation analysis,” in 52nd IEEE Conference on Decision and Control, Florence, 2013. 
[7]  B. Nicoletta, M. Davide and M. Riccardo, ”Topology Identification of Smart Microgrids,” 2013. 
[8]  M. H. Wen, R. Arghandeh, A. von Meier, K. Poolla and V. O. Li, ”Phase Identification in Distribution Networks with micro-synchrophasors,” in IEEE Power & Energy Society General Meeting, Denver, 2015. 
[9]  D. Deka, S. Backhaus and M. Chertkov, ”Estimating Distribution Grid Topologies: A Graphical Learning based Approach,” in Power Systems Computation Conference, Genova, 2016. 
[10]  R. Sonderegger, ”Smart grid topology estimator”. United States Patent US 2015 / 0241482 A1, 27 August 2015.
[11]  J. L. Kann, ”Grid topology mapping with voltage data”. United States Patent US 2016 / 0109491 A1, 21 April 2016.
[12]  X. Li, H. V. Poor and A. Scaglione, ”Blind Topology Identification for Power Systems,” in IEEE International Conference on Smart Grid Communications, Vancouver, 2013. 
[13]  V. Arya, T. S. Jayram, S. Pal and S. Kalyanaraman, ”Inferring connectivity model from meter measurements in distribution networks,” E-Energy, pp. 173-182, 2013. 
[14]  S. Jayadev, N. Bhatt, R. Pasumarthy and A. Rajeswaran, ”Identifying Topology of Power Distribution Networks Based on Smart Meter Data,” in https://arxiv.org/pdf/1609.02678, 2016. 
[15]  ”Grid diagram creation and opearation control”. United States Patent US 2015 / 0153153 A1, 4 June 2015.
[16]  V. Arya and R. Mitra, ”Voltage-based clustering to infer connectivity information in smart grids”. US Patent US20150052088A1, 19 February 2015.
[17]  J. Bickel and R. Carter, ”Automated hierarchy classification in utility monitoring systems”. US Patent US 20070005277 A1, 4 January 2007.
[18]  O. Xu, B. Geisser, A. Brenzikofer, S. Tomic and M. Eisenreich, ”Verfahren und einrichtung zur bestimmung der topologie eines stromversorgungsnetzes”. Europe Patent EP 3065250 A1, 7 September 2016.
[19]  V. Arya, S. Kalyanaraman, D. P. Seetharamakrishnan and J. Thathachar, ”Determining a connectivity model in smart grids”. US Patent US 20140039818 A1, 6 February 2014.
[20]  F. Kupzog and A. Lugmaier, ”Method for ascertaining parameters for power supply systems, and system for carrying out the method”. WO Patent WO 2013026675 A1, 28 February 2013.
[21]  ”Estimation Method of Distribution Network structure mixed integer quadratic programming Model”. China Patent CN 104600699 B, 17 August 2015.
[22]  ”Branch active power flow-based power network topology error identification method”. China Patent CN 201510973479 , 23 December 2015.
[23]  M. Coyne, N. Jespersen, P. Sharma, S. Srinivasan and L. Sheung Yin Tse, ”Historic storage of dual layer power grid connectivity model”. US Patent US 20120089384 A1, 12 April 2012.
[24]  Y. Sharon, A. Annaswamy, M. A. Legbedji and A. Chakraborty, ”Topology identification in distribution network with limited measurements”. US Patent US 20130035885 A1, 7 February 2013.
[25]  E. G. De Buda, ”System for accurately detecting electricity theft”. Canada Patent CA2689776, 12 July 2010.
[26]  P. Deschamps and M.-C. Alvarez-Herault, ”Method and device for determining the structure of a power grid”. Europe Patent EP2458340 A2, 30 May 2012.
[27]  V. Gouin, M.-C. Alvarez-Herault, P. Deschamps, S. Marie and Y. Lamoudi, ”Procede et systeme de determination de la structure d'un reseau de distribution d'electricite et programme d'ordinateur associe”. Europe Patent EP 3086093 A1, 26 October 2016.
[28]  V. Kekatos, G. B. Giannakis and R. Baldick, ”Online Energy Price Matrix Factorization for Power Grid Topology Tracking,” IEEE Transactions on Smart Grid , vol. 7, no. 3, pp. 1239 - 1248, 2016.
[29]  J. Kussyk, ”Low voltage distribution system and method for operating the same”, US Patent US 20140032144 A1, 30 Jan 2014.
【符号の説明】
【００６７】
１  電力ネットワーク
７  処理ユニット
１０  スラック分岐
Ｎ１  スラックノード
Ｎ２、Ｎ３  ノード
Ｌ１、Ｌ２  分岐
Ｍ１、．．．、Ｍ７  計測ユニット
【図１】
2020533926.tif 000054
【図２】
2020533926.tif 000055
【図３】
2020533926.tif 000056
【図４】
2020533926.tif 000057
【図５】
2020533926.tif 000058
【図６】
2020533926.tif 000059
【図７】
2020533926.tif 000060
【図８】
2020533926.tif 000061
【図９】
2020533926.tif 000062
【図１０】
2020533926.tif 000063
【国際調査報告】
2020533926.tif 000064
2020533926.tif 000065