JP2013506389A - 状態マトリックスに依存せずにｐｍｕ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法 - Google Patents

Abstract

本発明は、ネットワークのパラメータおよび状態マトリックスに依存せず、フェーザ計測ユニット（ＰＭＵ）を備えた測定点の実測動態データ、および電力系統データ収集と監視制御システム（ＳＣＡＤＡ）の測定データまたは状態推定データのみに基づいて、ＰＭＵを備えていない低可観測性測定点のダイナミクスをリアルタイムで推定する方法を開示する。当該方法は、再帰最小二乗法を利用して、一定期間における準定常状態のデータシーケンス、即ち連続した複数の断面のＳＣＡＤＡ測定データまたは状態推定データに基づいて、ＰＭＵを備えていないノードの推定しようとする量の変化量と、ＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合関係を見つける。そして、当該線形結合関係を利用して、ＰＭＵを備えていないノードの推定しようとする量の各ＰＭＵ測定時間断面における該当値を推定することにより、ＰＭＵを備えていない測定点に対するダイナミクスのリアルタイム推定を実現する。実際のデータによると、この方法は、非常に高い推定精度を実現することができ、工程応用上の誤差基準を満たしている。ＰＭＵの配置が不足している場合、ＰＭＵを備えていない測定点のダイナミクスを把握する必要がある時に、当該方法は大変実用的である。

Description

本発明は、電力系統の広域測定技術分野に属し、特に、ＰＭＵを備えていない測定点の測定量のダイナミクスを推定する方法に関する。

従来、電力系統コントロールセンターの情報は、主に、ＳＣＡＤＡシステム（データ収集および監視制御システム）から収集されたものであり、電圧、電流、電力、周波数などのリモート測定のサンプリング間隔が１〜５秒であり、しかも、タイムスケールがないため、電力ネットワークの準定常状態の情報しか反映することができない。近年来、グローバル・ポジショニング・システムであるＧＰＳに基づく同期フェーザ測定技術の電力系統への応用が益々広がり、当該技術に基づくコアデバイスは、各変電所および発電所に配置されるフェーザ測定ユニット、即ちＰＭＵである。ＰＭＵは、各変電所・発電所における、同期タイムスケールがマークされた電圧フェーザ、電流フェーザ、電力、周波数、発電機回転速度などの情報を、１０ｍｓ〜３０ｍｓのサンプリング間隔でコントロールセンターのマスターステーションに送信する。これによって、コントロールセンターにおいて、電力系統の動的挙動に対するオンライン集中監視が実現される。しかしながら、現在のＰＭＵ装置の価額がかなり高く、尚且つ、ＰＭＵのために大量な通信資源を装備する必要があるため、現段階および今後相当長期間にわたって、一部の変電所・発電所にしかＰＭＵを取り付けることができない。そのため、現在、電力ネットワークのコントロールセンターでは、ＰＭＵに基づく広域測定システム（ＷＡＭＳ）において、電力ネットワークの動的挙動に対する可観測性のすべての要求を満たすことがまだ実現されていない。

ＰＭＵを取り付けておらず（または電力会社に管轄されておらず）、しかも完全可観測性を有していない変電所・発電所またはノードのダイナミクスに対するリアルタイム監視を実現するために、補助オペレーターがオンラインで安全・安定対策を実施する。そのため、既存のオンライン測定情報を利用して、ＰＭＵを備えていないノード（以下、ＰＭＵ無しのノードと称することもある）のダイナミクスをリアルタイムで推定することは、非常に重要かつ実用的なことである。現在、様々な文献に開示された方法の中に、電力ネットワークにおける、ＰＭＵを備えたノード（以下、ＰＭＵ有りのノードと称することもある）の電圧フェーザとＰＭＵ無しのノードの電圧フェーザとの間の感度関係、および該当のＳＣＡＤＡ測定情報に基づいて、ＰＭＵ有りのノードの電圧フェーザの振幅値と位相角のダイナミクスによって、ＰＭＵ無しのノードの電圧フェーザの振幅値と位相角のダイナミクスを推定する方法が存在する。但し、ＰＭＵ有りのノードの電圧フェーザ（振幅値と位相角）とＰＭＵ無しのノードの電圧フェーザ（振幅値と位相角）との間の感度関係は、状態マトリックス、即ち潮流方程式のヤコビ行列によって導出される。このような方法では、以下の３つの問題点が存在する。
（１）原始ヤコビ行列の各要素が、ノード電圧の振幅値、位相角およびノードアドミタンスと関連しているため、潮流およびネットワーク状態が変化すると、それに伴ってヤコビ行列を常に更新する必要がある。しかしながら、ヤコビ行列更新の計算量が膨大であり、オンラインで１０ｍｓレベルのダイナミクスのデータを生成することは適切ではない。それに対して、電圧振幅値および位相角の変化を無視し、さらに、電力ネットワーク構成要素の抵抗を無視するなど一連の簡略化措置を取ることにより、ヤコビ行列の要素をノードのアドミタンス行列の虚部の要素に取り替えて、ＰＭＵ無しのノードとＰＭＵ有りのノードとの間の近似線形相関係数を求める方法が、ある文献に開示されている。なお、これらの相関係数は、ネットワークの構造が変わった時のみ更新され、その他の状況における潮流変化は、ノード間の相関係数の更新を引き起こさない。これによって、オンライン計算
の作業量は大幅に減少されるが、ダイナミクスのオンライン推定の精度も一定程度下がる。
（２）上記したように、簡略化されたヤコビ行列をアドミタンス行列から得ることはできるが、アドミタンス行列は、線路、変圧器、発電機などのネットワーク構成要素のアドミタンスパラメータに依存している。通常、これらのアドミタンスパラメータは、メーカが提供した定格値または典型的な値から得られたものである。これらの数値と実際のネットワーク構成要素のパラメータとの間は、時々、大きな誤差が存在する。しかも、これらの構成要素のパラメータは、温度、周波数などの運転環境または運転状態の変化に伴って変化することもある。これらのパラメータが精確でなければ、ＰＭＵ有りのノードとＰＭＵ無しのノードとの測定量の変化量間の線形相関関係を精確に計算することができず、さらに、ＰＭＵ無しのノードの測定量のダイナミクスを精確に推定することができない。また、電力ネットワークの構成要素が非常に多いので、上記パラメータが精確か否かについて、一つずつ検証してオンライン補正を行うことが不可能なため、ダイナミクスの推定が不精確になっても補正されない。
（３）既存の文献に開示されたダイナミクスの推定は、電圧フェーザのみに限られており、その他の測定、特に電力に対する処理方法に言及していない。

また、既存の文献において、ＰＭＵおよびＳＣＡＤＡデータに基づいて混合状態推定または混合動態状態推定を行うといった内容は多いが、いずれもＰＭＵデータを利用して、現在または将来のＳＣＡＤＡ測定断面の状態を推定する精度を向上するものであり、ＰＭＵ無しのノードの測定に対するダイナミクスの推定を行うという内容に関連していないと言える。

したがって、上記説明から分かるように、現在、ＰＭＵ無しの測定点の測定に対するダイナミクスの推定を行うための、精確で実用的な方法はまだ存在していない。このような状況において、本発明は、ネットワークのパラメータおよび状態マトリックスの更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、フェーザ測定ユニットであるＰＭＵを備えた測定点の実測動態データおよび電力系統データ収集と監視制御システムであるＳＣＡＤＡの測定データまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のダイナミクスをリアルタイムで推定する方法を提供する。本発明の方法は、線形化理論を理論基礎として、現在の断面とヒストリー断面のＳＣＡＤＡ（または状態推定）データに基づいて、最小二乗法を利用して、ＰＭＵ有りのノードの測定量の変化量とＰＭＵ無しのノードの対応測定量の変化量との間の線形結合関係を見つけ出し、そして、ＰＭＵ測定時間断面において、ＰＭＵ有りのノードの測定情報を用いて、ＰＭＵ無しのノードの該当する測定情報を推定する。ダイナミクス推定の全体のプロセスは、完全に実測データに依拠しており、電力ネットワークの構成要素のパラメータおよびヤコビ行列または状態マトリックスの更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がない。測定量の変化量間の線形関係が比較的強い場合、当該方法によるダイナミクス推定の誤差は非常に小さい。例えば、電圧フェーザのダイナミクスに対する推定の場合、通常、振幅値の誤差は±０．０５％以内、角度誤差は±０．０５°以内になっている。一方、測定量の変化量間の線形関係が比較的弱い場合、ダイナミクス推定の誤差は比較的大きいが、工程応用上の許容誤差範囲内に止まる。例えば、通常、低周波振動期間における電力ネットワークの有効電力に対するダイナミクス推定の誤差は、約±５％内に収まることが可能である。

ＰＭＵ無しの測定点のダイナミクスを推定する方法の精確性が、ネットワークのパラメータおよびヤコビ行列のリアルタイム更新に依存し、しかも、推定が電圧フェーザのみに限られているという従来技術の短所を克服するために、本発明は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、ＰＭＵ有りの測定点
の実測動態データおよびＳＣＡＤＡ測定データまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のダイナミクスをリアルタイムで推定する方法を開示する。当該方法では、ＰＭＵ無しの測定点は、完全可観測性を有する必要がない。本発明は、以下の具体的な技術案を用いる。

本発明に係るダイナミクスのリアルタイム推定方法は、一定期間における準定常状態のデータシーケンス、即ち、連続した複数の断面のＳＣＡＤＡデータまたは状態推定データに基づいて、ＰＭＵを備えていないノードの推定しようとする量の変化量と、ＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合関係を見つけ、そして、当該線形結合関係を利用して、ＰＭＵを備えていないノードの推定しようとする量のダイナミクスをリアルタイムで推定する。本発明を実現するための典型的なフローチャートは図１Ａと図１Ｂに示すようになっており、具体的な実施ステップは、以下の通りである。

第１のステップ：ダイナミクスの推定を行う必要があり、ＰＭＵを備えていないノードiを選択し、当該ノードにおける該当の構成要素の電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数および発電機回転速度などの、当該ノードにおけるダイナミクスを推定する必要がある物理量X_iを選択する。

第２のステップ：各ＰＭＵを備えていないノードi毎に、そのダイナミクスを推定するために使われ、ＰＭＵを備えたノードの集合S_PMU-iを選択し、当該ノードの集合は、通常、ＰＭＵを備えていないノードとの電気的距離が比較的近い複数の、ＰＭＵを備えたノードからなり、実用上、４〜６つのノードとすることが好ましい。本方法において、電気的距離に対する要求は厳しくなく、ＰＭＵを備えたノードと推定しようとするノードとが同じ電気島（electrical island）に属していることが保証される前提で、集合S_PMU-i中のノードを、オフラインで形成された典型的な運転方式におけるノードインピーダンスマトリックスに基づいて選択し、または手動で指定することができ、ＰＭＵを備えたノードとＰＭＵを備えていないノードとの間の最小線路段数に基づいて、オンライン自動探索方式で決定してもよい。具体的なステップとして、まず、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、ノードiと１段の線路を介して接続、即ち直接接続しているＰＭＵノードを集合S_PMU-iの中に優先に入れてから、集合中のノード数が予め決定された数になるまで、ノードiとそれぞれ２段の線路、３段の線路、４段の線路を介して接続しているＰＭＵを備えたノードを、介した線路の段数が小さい方から順に集合S_PMU-iに加える。集合S_PMU-i中のノードの、推定しようとする測定と同じ種類の測定をダイナミクス推定用に選択することができる。

第２のステップにおいて、次のデータ予備処理プロセスがさらに含まれている。

（ａ）推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵを備えたノードの電圧振幅値または電流振幅値を、ＰＭＵを備えていないノードの電圧レベルにおける電圧振幅値または電流振幅値に変換しなければならない。

（ｂ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するあるノードの電圧フェーザを基準とした相対位相角に変換し、当該相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算するようにしなければならない。

（ｃ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、変圧器の
結線方式によって生じた位相シフトに基づいて、ＰＭＵを備えたノードの電圧位相角または電流位相角を、ＰＭＵを備えていないノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算するようにしなければならない。

第３のステップ：各ＰＭＵを備えていないノードiにおける推定しようとする測定量X_iの変化量ΔX_iと該当のＰＭＵを備えたノードにおける対応の測定量X_jの変化量ΔX_jとの間の線形結合係数b_ij、および再帰最小二乗法パラメータに対して、初期化を行い、但し、j∈S_PMU-iとする。初期化の具体的な方法は以下の通りである。

（ａ）まず、推定しようとする測定量の変化量ΔX_iの線形結合係数b_ijを求めるという問題の再帰最小二乗法の数学モデルを作成する。再帰最小二乗法を用いて問題を解く場合の一般的な形式は、

であり、但し、dは観測された出力ベクトル、wは求めようとするパラメータベクトル、Aは複数のベクトルu_jからなる状態マトリックスであり、u_jはパラメータベクトルwと観測量d( j )との関係を表すものであり、vはノイズベクトルである。推定しようとする測定量の変化量ΔX_iの線形結合係数b_ijを求めるという問題について、上式における計算中に必要な変数の定義は、

となっており、ΔX_iは、推定しようとする測定量X_iの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との間の差値であり、各推定しようとする測定をそれぞれ推定するため、実際には、出力ベクトルdは１つの要素しかない。ΔX₁ ΔX₂ …ΔX_j …ΔX_N は、それぞれ、S_PMU-i中の各ノードの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値であり、b_i1 b_i2 … b_ij … b_iNi は、それぞれ、S_PMU-i中の各ノードの該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔX_iとの間の線形結合係数である。

（ｂ）前記推定しようとする測定量の変化量ΔX_iの線形結合係数 b_ijの初期値を、通常1/N_iとすると、再帰最小二乗法で求めようとするパラメータベクトルwの初期値はw_-1=[1/N_i 1/N_i … 1/N_i ]^Tとなり、但し、N_iはＰＭＵを備えていないノードiのS_PMU-i集合中の
要素の数であり、下付き「-1」は初期値を表す。

（ｃ）前記各線形結合係数と対応する、ＰＭＵを備えた測定の変化量に対して、重み付け係数を設定することにより、N_i×N_i次の対角マトリックスである重みマトリックスΠが得られ、その対角要素は、各重み付け係数からなる。ＰＭＵを備えた測定の変化量の測定品質が高いほど、また、そのノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きい。再帰最小二乗法の計算プロセスにおける中間プロセスマトリックス変数Pの初期値P_-1をΠとし、ＰＭＵを備えた測定の変化量の測定品質が不明な場合、各重み付け係数を同じ値の0.1とし、この場合、

となり、但し、下付き「-1」は初期値を表す。

（ｄ）望ましくは、初期化忘却因子λを０．９９とする。

第４のステップ：各推定しようとする測定量X_i毎に、現在最新のＳＣＡＤＡ測定時刻T_kにおける計測値または最新の状態推定時刻T_kにおける計算値、およびヒストリー時刻におけるＳＣＡＤＡ計測値または状態推定値に基づき、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔX_i(T_k)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔX_j(T_k)との間の線形結合係数b_ij(T_k)を算出し、但し、j∈S_PMU-iとし、T_kの下付き「k」は第k回のＳＣＡＤＡ測定または第k回の状態推定を示し、循環変数kは0から始まる整数を示す。具体的には、以下のステップ（ａ）〜（ｃ）を含む。

（ａ）ノードiの推定しようとする測定の現在の第k時間ステップ、即ちT_k時刻におけるＳＣＡＤＡ測定値または状態推定値の、前回のサンプリング時刻におけるＳＣＡＤＡ測定値または状態推定値である第k-1時間ステップにおける対応値に対する変化量を算出することにより、下記第k時間ステップの出力ベクトルd_kが得られる。

（ｂ）次の式により、集合S_PMU-i中の各ノードの、推定しようとする測定と対応する各測定の現在の第k時間ステップにおけるＳＣＡＤＡ測定値または状態推定値の、第k-1時間ステップにおける対応値に対する変化量のベクトルu_kを算出する。

（ｃ）次の３つの回帰公式に基づいて、第k時間ステップにおける再帰最小二乘推定法
の中間プロセスマトリックスg_k、パラメータベクトルマトリックスw_k=[ b_i1(T_k) b_i1(T_k) … b_iNi(T_k) ]^T、中間プロセスマトリックスP_kを順に算出する。循環変数kは第0時間ステップから始まり、第0時間ステップの計算に使われたP_-1およびw_-1は、第３のステップにおける初期化処理プロセスによって得られるものである。第k時間ステップにおける各線形組み合わせ係数b_ij(T_k)は、第k時間ステップで得られたw_kによって得られる。

第５のステップ：第４のステップにおいて、T_k時刻におけるＳＣＡＤＡ計測値または状態推定値に基づいて算出された、現在推定しようとする量の変化量ΔX_iとＰＭＵを備えたノードにおける対応の測定量の変化量ΔX_jとの間の線形結合係数b_ij(T_k)に基づいて、以下の式

を用いて、次回のＳＣＡＤＡ測定時刻T_k+1または状態推定時刻T_k+1までの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵを備えていないノードの推定しようとする量X_i(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻T_kにおける該当の量X_i(T_k)に対する変化量の推定値ΔX_i(t)を算出し、但し、T_k<t<T_k+1とする。

さらに、方程式X_i(t)=ΔX_i(t)+X_i(T_k) に基づいて、当該ＰＭＵを備えていないノードの推定しようとする量の各ＰＭＵ測定時刻における推定値X_i(t)を算出し、但し、T_k<t<T_k+1とする。なお、推定された相対位相角の値を-180°〜+180°の範囲内に換算し、または、推定された相対位相角の値を絶対位相角の値に変換してから、-180°〜+180°の範囲内に換算するようにしなければならない。

ＳＣＡＤＡ測定または状態推定値が存在する時刻において、ＰＭＵ無しのノートの推定しようとする量は、ＳＣＡＤＡ測定値または状態推定値とする。

さらに、第５のステップにおいて、次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間に、電力ネットワークに短絡または断線事故が発生する可能性があり、事故が発生した場合、ＳＣＡＤＡの該当スイッチング量が即時に更新されるが、アナログ量は、通常、即時に更新されない。この場合、ネットワーク構成の変化により、ヤコビ行列またはノードアドミタンス行列が変わるため、ＰＭＵを備えた各ノードによる、ＰＭＵを備えていないノードの推定量の変化への影響は大きく変化する。しかしながら、ＳＣＡＤＡのアナログ測定データまたは状態推定データの更新が遅いため、新しい係数を直ぐに生成することが難しい。このような推定上の誤差を減らすために、トポロジー解析によって、ＰＭＵを備えていないノードiまたはＰＭＵを備えたノードj（j∈S_PMU-i）に直接接続されている線路に、短路または断線が発生した場合、以下の処理（ａ）〜（ｃ）を実施する。即ち、
（ａ）ＰＭＵを備えたノードjが系統から孤立された場合、それと対応した線形結合係
数b_ijをその他のＰＭＵを備えたノードに平均配分し、孤立されたＰＭＵを備えたノードjの線形結合係数b_ijを0とし、
（ｂ）ＰＭＵを備えていないノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において当該ノードのダイナミクスを推定することが不可能であることを示し、孤立されたＰＭＵを備えていないノードiが再び系統に接続されたら、孤立される前の線形結合係数を初期係数b_ijとして、該当の測定のダイナミクスを推定し、
（ｃ）その他の状況、即ち、短絡または断線によってＰＭＵを備えていないノードiまたはＰＭＵを備えたノードjが系統から孤立されたという状況にならなかった場合、特殊な処理を行わず、次回の、即ち第k+1時間ステップのＳＣＡＤＡアナログ測定データまたは状態推定データの更新によって、線形結合係数b_ijを補正して、ネットワーク構成の変化を表す。

第６のステップ：新しいＳＣＡＤＡ測定データまたは状態推定結果が得られると、第４のステップに戻り、ＰＭＵを備えていないノードの推定しようとする量の変化量とＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合係数を更新する。

本発明は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、ＰＭＵ有りの測定点の実測動態データおよびＳＣＡＤＡデータまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のＰＭＵ測定期間におけるダイナミクスのリアルタイム推定を実現することができ、従来のダイナミクス推定方法における、精確性が電力ネットワークパラメータおよびヤコビ行列の精確性に依存しているため、計算量が大きく、リアルタイム性が低く、精確性が確保できないといった問題点を克服することができる。本発明の方法を利用する場合、測定量の変化量間の線形関係が比較的強い時に、ダイナミクス推定の誤差は非常に小さい。例えば、電圧フェーザに対するダイナミクス推定の場合、通常、その振幅値の誤差は±０．０５％以内、角度の誤差は±０．０５°以内である。一方、測定量の変化量間の線形関係が比較的弱い時に、ダイナミクス推定の誤差は相対的に大きいが、工程応用上の許容可能な誤差範囲内に止まる。例えば、低周波振動期間における電力ネットワークの有効電力に対するダイナミクスの推定の場合、通常、その誤差を約±５％以内に保つことができる。なお、ＳＣＡＤＡまたは状態推定データ断面間の時間間隔の短縮に伴って、上述したダイナミクス推定の誤差がさらに小さくなる。

状態マトリックスに依存せずにＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法のフローチャートである。 状態マトリックスに依存せずにＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法のフローチャートである。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（近くのＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が１秒である）。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（近くのＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が５秒である）。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（離れたＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が１秒である）。 ５００ＫＶ側測定によって光谷ステーション４号変圧器２２０ＫＶ母線の電圧振幅値ダイナミクスを推定した場合の相対誤差グラフである。 光谷ステーション夙光線の電圧位相角ダイナミクスの推定値と実測値との絶対誤差グラフである。 光谷ステーション夙光線の有効電力ダイナミクスの推定値と実測値との相対誤差グラフである。 ある低周波振動の前および振動中に、非協同領域の振動電力を用いて万龍１線の有効電力ダイナミクスを推定した場合の相対誤差グラフである。

本発明に開示される方法は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、ＰＭＵ有りの測定点の実測動態データおよびＳＣＡＤＡ測定データまたは状態推定データのみに基づいて、ＰＭＵ無しのノードの測定量のダイナミクスに対するリアルタイム推定を実現することができる。本発明は、ある測定量の変化量を、当該測定量と相関する複数の測定量の変化量の線形結合から得ることができるという発想に基づくものである。即ち、

上式において、Xiはノードiにおける、ダイナミクスが推定される１つの測定量であり、電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数または発電機回転速度などであってもよい。ΔX_i(t)は、時刻tにおけるX_iの変化量を表すものである。集合S_PMU-iは、例えば、変電所および発電所のトポロジーノードのような、ノードiとの電気的距離が近いＰＭＵ有りのノードからなる。S_PMU-i中のノードに対して、ノードiとの電気的距離が近いことを要求しているが、この要求は厳格な要求ではない。後述の計算例から分かるように、ノードiとの電気的距離が比較的遠いが、ノードiと同一電気島内に存在するノードを利用しても、同様に、ノードiのダイナミクスに対する推定を実現することができる。この場合、誤差は大きくなるが、通常の工程応用上の誤差許容範囲内に収まる。X_jは集合S_PMU-i中のノードjにおける、X_iと同じ性質の測定量であり、b_ijは相関する測定量の変化量の線形結合係数である。

本方法において、集合S_PMU-iの大きさに関する厳密な規定はない。式（１）から分かるように、S_PMU-i中の要素が1になってはいけない。さもなければ、推定されたノードのダイナミクスの形状は、ＰＭＵ有りのノードのダイナミクスと完全に一致してしまう。様々な要因による、推定されるノードのダイナミクスの形状および大きさへの影響を反映するために、集合S_PMU-i中の要素として、数が多く、分布範囲が広いことが望ましい。しかしながら、集合S_PMU-i中の要素の増加に伴って、ダイナミクスのリアルタイム推定の計算量も増加するため、ダイナミクス推定のリアルタイム性に影響を与える可能性がある。実際の工程応用において、集合中の要素の数を４〜６とすることが最も適切であると実験によって判明されている。

また、式（１）に示した線形結合関係は、主に経験に基づいた公式であり、厳密な数学関係に依存するものではないため、測定量の変化量間の線形関係が著しい電圧フェーザのダイナミクス推定に応用することができ、電力、周波数など変化量間の線形関係が著しくない測定量のダイナミクス推定に応用することもできる。ユーザは、線形化する具体的な数学モデルに関わる必要がなく、再帰最小二乗法を用いてヒストリーデータ学習を行うことによって、線形化係数を求めることができる。実験で証明されたように、本発明の方法
によれば、ダイナミクス推定に関与する測定量の変化量の間に厳密な線形関係が存在する場合、動態推定結果の精度は非常に高く、例えば、電圧振幅値のダイナミクス推定の誤差は、通常、±０．０５％以内であり、電圧位相角度のダイナミクス推定の誤差は、通常、±０．０５°以内である。一方、ダイナミクス推定に関与する測定量の変化量間の線形関係が比較的に弱い場合、ダイナミクス推定の誤差は比較的に大きいが、工程応用上の誤差許容範囲内に収まり、例えば、低周波振動期間における電力ネットワークの有効電力のダイナミクスを推定する場合でも、通常、誤差を約±５％内に維持することができる。なお、ＳＣＡＤＡまたは状態推定データ断面間の時間間隔の短縮に伴い、例えば、５秒から１秒に変更すると、上記ダイナミクス推定の誤差はさらに小さくなる。

式（１）および上述した原理に基づいて、連続した複数の断面のＳＣＡＤＡ測定データまたは状態推定データを利用して、ＰＭＵ無しのノードの推定しようとする測定量の変化量とＰＭＵ有りのノードにおける対応の測定量の変化量との間の線形結合関係を求める。そして、当該線形結合関係を利用して、ＰＭＵ無しのノードの推定しようとする測定量の各ＰＭＵ測定時間断面における該当値を推定することにより、ＰＭＵ無しのノードの測定量のダイナミクスに対するリアルタイム推定を実現する。

準定常状態のデータシーケンスのソースに関して、電圧振幅値、電流振幅値、有効電力、無効電力および周波数などの電気量のみに対してダイナミクスを推定する場合は、ＳＣＡＤＡデータのみを使用すればよく、通常、そのデータ断面間の間隔は１〜５秒であり、さらに、電圧位相角、電流位相角のダイナミクスを推定する必要がある場合では、状態推定の結果を使用しなければならない。その理由は、現在、ＳＣＡＤＡが位相角の測定を行っていないからである。状態推定データ断面間の時間間隔は、状態推定プログラムの計算効率および電力ネットワークの規模によって異なり、ダイナミクスのリアルタイム推定の応用において、通常、５秒に設定することができる。また、現在、一部のＳＣＡＤＡおよびＰＭＵ／ＷＡＭＳ（広域測定システム）は、発電機の回転速度など非電気量の収集もサポートしている。この場合では、本発明の方法を用いて、ＰＭＵ無しの測定点における非電気量に対するダイナミクス推定を実現することができる。

式（１）に対して、再帰最小二乗法を用いて係数b_ijを求めるようにするには、次の理由がある。第１、再帰の方法によれば、新しいＳＣＡＤＡまたは状態推定断面を得た直後に当該時間断面の係数b_ijを求めることができるため、計算のリアルタイム性を維持することができる。第２、式（１）は、以下に示す、再帰最小二乗法で問題を解く場合の一般的な形式に適合している。

式（２）において、dは観測された出力ベクトル、wは求めようとするパラメータベクトル、Aは複数のベクトルu_jからなる状態マトリックスであり、u_jはパラメータベクトルwと観測量d( j )との関係を表すものであり、vはノイズベクトルである。したがって、式（１）と式（２）とは、以下の対応関係を有している。

但し、ΔX_iは、推定しようとする測定量X_iの現在のＳＣＡＤＡ測定値（または状態推定値）と前回のＳＣＡＤＡ測定値（または状態推定値）との差値を表すものであり、各推定しようとする測定をそれぞれ推定するため、実際に、出力ベクトルdは1つの要素しかない
。ΔX₁ ΔX₂ … ΔX_j … ΔX_Ni は、それぞれ、S_PMU-i中の各ノードにおける現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値である。

b_i1 b_i2 … b_ij … b_iNi は、それぞれ、S_PMU-i中の各ノードにおける該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔX_iとの間の線形結合係数である。式（２）は再帰最小二乗法で問題を解く場合の一般的な形式に適合している。そのため、再帰最小二乗法による汎用解法によって、係数b_ijを便利に算出することができる。

実際のシステムにおいて、本発明の原理に基づくダイナミクスリアルタイム推定ソフトウェアは、以下のような典型的なハードウェア環境で実行することが可能である。コントロールセンターの広域測定システム（ＷＡＭＳ）のハイレベル応用サーバまたはリアルタイムデータサーバにおいて、ダイナミクスリアルタイム推定ソフトウェアを実行することができる。ダイナミクスリアルタイム推定ソフトウェアは、ＳＣＡＤＡ／ＥＭＳシステムのリアルタイムデータベースから、現在の電力ネットワークのＳＣＡＤＡ測定値または状態推定値をオンラインで獲得し（推定しようとする量およびそれと相関するＰＭＵ有りのノードにおける対応の測定量のＳＣＡＤＡ測定値または状態推定値のみを獲得してもよい）、そして、再帰最小二乗法を利用して、現在のＳＣＡＤＡまたは状態推定時間断面に対応する各推定しようとする測定量の各線形結合係数を計算する。次回のＳＣＡＤＡ測定値または状態推定値を受信するまで、ダイナミクスリアルタイム推定ソフトウェアは、現在の線形結合係数およびＷＡＭＳリアルタイムデータサーバから獲得した関連のＰＭＵ有りのノードにおける対応の動態測定値を利用して、ＰＭＵ無しのノードの推定しようとする測定量の各ＰＭＵ測定時間断面における該当値を推定し、そして、各ＰＭＵ測定時間断面における推定結果をＷＡＭＳリアルタイムデータサーバに保存し、さらに、ＰＭＵ無しの実測値であることを示すマークを付ける。

次に、実際のデータを用いて、本発明の技術案の具体的なステップをさらに説明する。図１Ａ、図１Ｂは、本発明の方法をソフトウェアで実現するためのフローチャートを示している。本発明の方法の有効性を検証するために、以下の計算例において、実際にＰＭＵで測定することができる測定量を、ダイナミクスを推定しようとする測定量とし、そして、推定したダイナミクスと実際に測定したダイナミクスとを比較することにより、本発明の方法の有効性を検証するようにしている。ネットワークパラメータおよび状態マトリックスに依存せずに、ＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定するアルゴリズムは、以下のステップによって実現される。

（１）ダイナミクスを推定する必要があるＰＭＵ無しのノード、および当該ノードにおける、ダイナミクスを推定する必要がある物理量X_iを選択し、当該物理量X_iは、例えば、ある構成要素の電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効、無効、周波数または発電機回転速度などである。

（２）各ＰＭＵ無しのノードi毎に、そのダイナミクスを推定するために使われる、ＰＭＵ有りのノードの集合S_PMU-iを選択する。当該ノードの集合は、通常、ＰＭＵ無しのノードとの電気的距離が比較的近い複数の、ＰＭＵ有りのノードからなり、手動指定の方法または自動探索の方法によって決定される。自動探索方法の１つとして、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、１段の線路を介してノードiと接続、即ち直接接続しているＰＭＵノードを集合S_PMU-iに優先に入れて、その後、集合中のノードの数が４〜６（状況によって集合の大きさを選ぶ）になるまで、それぞれ、２段の線路、３段の線路、４段の線路などを介してノードiと接続しているＰＭＵ有りのノードを、介した線路の段数が小さいほうから順に集合S_PMU-iに加える。集合S_PMU-i中のノードの、推定しようとする測定量と同じ種類の測定をダイナミクスの推定用に選択することができる。推定される量によって、原始測定データに対して、以下の予備処理を行う必要がある。

（２．ａ） 推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵ有りのノードとＰＭＵ無しのノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵ有りのノードの電圧振幅値または電流振幅値を、ＰＭＵ無しのノードの電圧レベルにおける電圧振幅値または電流振幅値に変換しなければならない。

（２．ｂ） 推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するあるノードの電圧フェーザを基準とした相対位相角に変換し、即ち、相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算するようにしなければならない。

（２．ｃ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵ有りのノードとＰＭＵ無しのノードとは電圧レベルが異なっているならば、変圧器の結線方式で生じた位相シフトに基づいて、ＰＭＵ有りのノードの電圧位相角または電流位相角を、ＰＭＵ無しのノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算するようにしなければならない。

（３）各ＰＭＵ無しのノードiにおける推定しようとする測定の係数b_ijおよび再帰最小二乗法のパラメータに対して、初期化を行う。

（３．ａ） 各推定しようとする測定の係数b_ijの初期値は、1/N_iとすることができる。N_iはＰＭＵ無しのノードiのS_PMU-i集合中の要素の数である。同様に、再帰最小二乗法の未知数ベクトルwに対する初期化を行って、w_-1=[1/N_i 1/N_i … 1/N_i ]^Tが得られ、下付き「-1」は初期値を表すものである。

（３．ｂ） 各線形結合係数に対応したＰＭＵ有りの測定の変化量に重み付け係数を設定して、N_i×N_i次の対角マトリックスである重みマトリックスΠが得られる。その対角要素は、各重み付け係数からなる。ＰＭＵ有りの測定の変化量の測定品質が高いほど、また、その属するノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きい。再帰最小二乗法の計算プロセスにおける中間プロセスマトリックスの変数Pの初期値P_-1はΠとし、その下付き「-1」は初期値を表すものである。ＰＭＵ有りの測定の変化量の測定品質が判明されていない場合、各重み付け係数をいずれも同じ値の０．１とする。この場合、

（３．ｃ） 初期化忘却因子λは、０．９９とすることができる。

（４）各推定しようとする測定量X_i毎に、以下の繰り返し処理を実行する。各最新時刻（第k時間ステップ）の準定常状態断面において、対応のＳＣＡＤＡ測定値または状態推定値を利用して、以下のステップに従って係数b_ijを計算し、さらに、ＰＭＵ有りのノードの動態測定データを利用して、次回の準定常状態時間断面までの期間中の推定される量のダイナミクス（グラフ）を推定する。

（４．ａ）ノードiの推定しようとする測定の現在の第k時間ステップ（即ち、T_k時刻）におけるＳＣＡＤＡ計測値または状態推定値の、前回のＳＣＡＤＡまたは状態推定断面である第k-1時間ステップにおける対応値に対する変化量を算出することにより、第k時間ステップにおける出力ベクトルd_kを獲得する。即ち、

（４．ｂ）変化量ベクトルu_kを算出する。但し、変化量ベクトルu_kは、集合S_PMU-i中の各ノードの、推定しようとする測定と対応する各測定の現在の第k時間ステップにおけるＳＣＡＤＡ計測値または状態推定値の、第k-1時間ステップにおける対応値に対する変化量のベクトルu_kである。即ち、

（４．ｃ）次の式に従って、マトリックスg_kを求める。

（４．ｄ）次の式に従って、ベクトルw_k=[ b_i1(T_k) b_i1(T_k) … b_iNi(T_k) ]^Tを求める。即ち、各線形結合係数b_ij(T_k)を求める。

（４．ｅ）次の式に従って、マトリックスP_kを求め、次のＳＣＡＤＡまたは状態推定時刻T_k+1における線形結合係数計算のための準備を行う。

i） 各ＰＭＵ測定時刻tにおいて、推定しようとする量X_iのS_PMU-i集合中の各ノードのＰＭＵ測定量の、ＳＣＡＤＡまたは状態推定の第k時間ステップの値に対する変化量であるΔX_j(t)=X_j(t)−X_j(T_k) を算出し、但し、T_k<t<T_k+1。

ii）ステップ（４．ｄ）で算出したＳＣＡＤＡまたは状態推定の第k時間ステップの線形結合係数b_ij(T_k)を用いて、次の式

に基づいて、ＳＣＡＤＡまたは状態推定の第k+1時間ステップになるまでの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵ無しのノードの推定しようとする量X_i(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻T_kおける該当量X_i (T_k)に対する変化量の推定値ΔX_i(t)を算出する。但し、T_k<t<T_k+1 である。

iii）方程式X_i(t)=ΔX_i(t)+X_i(T_k) に基づいて、当該ＰＭＵ無しのノードの推定しようとする量の、各ＰＭＵ測定時刻における推定値X_i(t)を算出し、但し、T_k<t<T_k+1。推定された相対位相角値を-180°〜+180°の範囲内に換算しなければならない。必要に応じて、推定された相対位相角値をさらに絶対位相角値に変換することができる。同様に、絶対位相角値も-180°〜+180°の範囲内に換算しなければならない。

ＳＣＡＤＡ測定または状態推定値が存在する時刻において、ＰＭＵ無しのノードの推定しようとする量は、ＳＣＡＤＡ測定値または状態推定値とする。

iv）現在のＰＭＵ測定時刻における推定しようとする量の推定値X_i(t)を記憶し、動態グラフを更新する。

v）現在のＰＭＵ測定時刻において、その他の推定しようとする量があるか否かを判断し、あると判断した場合は、ステップi）に戻り、ないと判断した場合、引き続きステップvi）を実行する。

vi）次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間、即T_k+1時刻になるまでの期間中、ＳＣＡＤＡのスイッチング量が、電力ネットワークの短路または断線故障のため、即時に更新されていることがあったか否かを検出する（通常、ＳＣＡＤＡのアナログ量は、短路または断線事故の場合に即時に更新されることがない）。上記更新があった場合、トポロジー解析を行い、ＰＭＵ無しのノードiまたはＰＭＵ有りのノードj（j∈S_PMU-i）と直接接続されている線路に短路または開路が発生した場合、以下の処理を行う。

（ａ）ＰＭＵ有りのノードjが系統から孤立された場合、それと対応した線形結合係数b_ijをその他のＰＭＵ有りのノードに均等配分し、孤立されたＰＭＵノードjの線形結合係数を０にする。

（ｂ）ＰＭＵ無しのノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において、当該ノードのダイナミクスを推定することができないことを示す。孤立されたＰＭＵ無しのノードが再び系統に接続されたら、孤立される前の線形結合係数を初期係数b_ijとして、該当の測定のダイナミクスを推定する。

（ｃ）その他の状況において、即ち、ＰＭＵ無しのノードiまたはＰＭＵ有りのノードj（j∈S_PMU-i）が、短路または断線のため、系統から孤立されたという事態が発生しなかった場合、特殊な処理を行わない。次回、即ち第k+1時間ステップのＳＣＡＤＡアナログ量測定または状態推定データの更新のみによって線形結合係数b_ijを修正して、ネットワーク構成の変化を反映する。

vii）新しい断面、即ち第k+1時間ステップのＳＣＡＤＡまたは状態推定断面が現れたか否かを判断し、現われていない場合、ステップi）に戻って、引き続き、各推定しようとする量X_iの既存の線形結合係数b_ij(T_k)を用いて、次のＰＭＵ測定時刻tにおける各推定しようとする量の数値を算出し、一方、現われた場合では、ステップ（４．ｆ）の繰り返しを中断して、ステップ（４．ａ）に戻って、各推定しようとする量X_iの新しい断面、即ち第k+1 時間ステップのＳＣＡＤＡまたは状態推定断面の線形結合係数b_ij(T_k+1)の計算を開始する。

次に、複数の実例を用いて、本発明の方法の有効性を検証する。本発明によるダイナミクスの精度を検証するため、選択した推定しようとする量は、実際にＰＭＵ測定が可能な量としている。そのため、推定されたダイナミクスを実際のＰＭＵ測定によって得られたダイナミクスと比較することにより、相対誤差グラフまたは絶対誤差グラフ（位相角測定精度を測るために使われる）を算出して、本発明の有效性を測ることができる。但し、ＰＭＵデータのサンプリング間隔を２０msec、ＳＣＡＤＡデータのサンプリング間隔を１secまたは５secとしている。本願の明細書において、相対誤差の計算公式は
相対誤差=（推定値−実測値）／実測値 （３）
であり、絶対誤差の計算公式は
絶対誤差=推定値−実測値 （４）
である。

実例１：ある省の光谷ステーション付近のＰＭＵ有りの変電所の電圧を用いて、光谷ステーションの５００ＫＶ側の夙光線の前端電圧のダイナミクスをリアルタイムで推定する。採用されたＰＭＵ有りの変電所は、慈湖ステーション、大吉ステーション、咸寧ステーション、玉賢ステーションである。これらの４つのステーションにおける、５００ＫＶ側の該当線路の前端電圧のＰＭＵ測定を選んで、ダイナミクスの推定を行う。これらのステーションにおける測定ノードは、光谷ステーション夙光線の前端ノードと直接接続され、または１つの他の変電所ノードのみを介して接続されている。ＳＣＡＤＡ測定のサンプリング間隔が１秒である場合、再帰最小二乗法を利用して推定した最初の３００秒の電圧振幅値のダイナミクスの推定値と実測値との相対誤差グラフは、図２に示すようになっている。相対誤差グラフから分かるように、ダイナミクス推定の相対誤差範囲は（-0.0456%、0.0399%）である。

実例２：上述した例と同じ測定環境において、ＳＣＡＤＡ断面の時間間隔を５秒に変更した場合、得られた夙光線の電圧振幅値の推定値と実測値との相対誤差グラフは、図３に示すようになっている。相対誤差グラフから分かるように、ダイナミクス推定の相対誤差範囲は（-0.11%、0.14%）である。誤差範囲は、１秒のサンプリング周期の場合に比べて、２〜３倍拡大されているが、最初からの６０秒、即ち１２つのＳＣＡＤＡサンプリングポイントの後に、誤差がほぼ（-0.05%、0.05%）以内に低減されている。再帰最小二乗法をオンラインで使用する場合、連続で使用することが考えられるため、計算の初期段階の誤差の影響を無視することができる。したがって、実際に応用する場合、誤差はほぼ（-0.05%、0.05%）以内に収まる。

実例３：光谷ステーションから比較的離れたＰＭＵ有りの変電所の電圧を用いて、光谷ステーション５００ＫＶ側の夙光線の前端電圧のダイナミクスをリアルタイムで推定する。採用されたＰＭＵ有りの変電所は、大別山ステーション、双河ステーション、孝感ステーション、興隆ステーションである。これらの４つのステーションにおける５００ＫＶ側の該当線路の前端電圧のＰＭＵ測定を選んで、ダイナミクスの推定を行う。これらのステーションにおける測定ノードと、光谷ステーションの夙光線の前端ノードとの間に、１〜４つの他の変電所ノードが介在している。ＳＣＡＤＡのサンプリング間隔が１秒である場合、再帰最小二乗法を利用して推定した最初の３００秒の電圧振幅値のダイナミクスの推定相対誤差グラフは、図４に示すようになっており、誤差範囲が（-0.0520%、0.0639%）である。

実例４：５００ＫＶのＰＭＵ測定を用いて、光谷ステーションの４号変圧器の２２０ＫＶ側の電圧のダイナミクスをリアルタイムで推定する。採用されたＰＭＵ有りの変電所は、慈湖ステーション、大吉ステーション、咸寧ステーション、玉賢ステーションである。これらの４つのステーションにおける５００ＫＶ側の該当線路の前端電圧のＰＭＵ測定を選んで、ダイナミクスの推定を行う。これらのステーションにおける測定ノードと、光谷ステーション４号変圧器の２２０ＫＶ側母線との間に、１〜４つの他の変電所ノードが介在している。ＳＣＡＤＡ周期が１秒である場合のダイナミクス推定の相対誤差のグラフは、図５に示すように、計算開始から２０秒を経過した後、相対誤差の範囲は（-0.10%、0.0997%）、即±0.1%になっている。再帰最小二乗法をオンラインで使用する場合、連続で使用することが考えられるので、計算の初期段階の誤差の影響を無視することができる。したがって、実際に応用する場合の誤差が（-0.1%、+0.1%）以内であると思われる。

実例５：ある省の光谷ステーション付近のＰＭＵ有りの変電所の電圧位相角を用いて、光谷ステーションの５００ＫＶ側の夙光線の前端電圧位相角のダイナミクスをリアルタイムで推定する。採用されたＰＭＵ有りの変電所は、慈湖ステーション、大吉ステーション、咸寧ステーション、玉賢ステーションである。これらの４つのステーションにおける５００ＫＶ側の該当線路の前端電圧位相角のＰＭＵ測定を用いて、ダイナミクスの推定を行
う。ＳＣＡＤＡのサンプリング間隔が１秒である場合、再帰最小二乗法を利用して推定された最初の３００秒の電圧位相角のダイナミクスの推定値と実測値との絶対誤差は、図６に示すようになっている。絶対誤差グラフから分かるように、当該電圧位相角に対するダイナミクス推定の絶対誤差範囲は（-0.043°、0.0426°）である。

実例６：ある省の光谷ステーション付近のＰＭＵ有りの変電所の線路の有効電力を用いて、光谷ステーション５００ＫＶ側の夙光線前端の線路の有効電力のダイナミクスをリアルタイムで推定する。採用されたＰＭＵ有りの変電所は、慈湖ステーション、大吉ステーション、咸寧ステーション、玉賢ステーションである。これらの４つのステーションにおける５００ＫＶ側の該当線路前端の有効電力のＰＭＵ測定を用いて、ダイナミクスの推定を行う。ＳＣＡＤＡサンプリング間隔が１秒である場合、再帰最小二乗法を利用して推定された最初の３００秒の線路有効電力のダイナミクスの推定値と実測値との相対誤差は、図７に示すようになっている。相対誤差グラフから分かるように、当該有効電力に対するダイナミクス推定の相対誤差は（-1.95%、2.72%）以内である。

実例７：実際に発生した低周波振動のケースを使って、本発明の方法による電力のダイナミクスのリアルタイム推定の能力を検証する。当該低周波振動の発生前および発生期間中の５００ＫＶ万龍１号線の有効電力のダイナミクスを推定する際に、極端に不利な状況における本発明の方法の推定精度を検証するために、特別に、万龍１号線と同じ協同領域に存在していない４つのＰＭＵ測定有りの５００ＫＶ線路のＰＭＵ実測有効電力を選んで、万龍１号線の有効電力のダイナミクスを推定した。得られた相対誤差グラフは図８に示すようになっている。なお、選んだ４つの線路は、長南１号線、葛岡線、江復１号線、江陵１号線である。電力ネットワークの外乱は２５６．５８秒で発生した。その後、０．４７Ｈｚを主モードとする電力ネットワークの低周波振動が引き起こされ、図８から分かるように、外乱発生後の極めて短い期間（０．１秒未満）内および再帰最小二乗法が起動した直後の一部の時刻での有効電力のダイナミクス推定の誤差は比較的大きいが、それ以外の期間において、低周波振動があるか否かにかかわらず、有効電力のダイナミクス推定の誤差は、すべて（-5%、+5%）の範囲内に収まっている。

上述した実例から、本発明に関する以下の検証結果が得られる。

（１）実例１および実例５から分かるように、測定量の変化量間の線形関係が比較的強い場合、本発明の方法を用いたダイナミクスのリアルタイム推定の誤差は非常に小さい。例えば、電圧フェーザのダイナミクスを推定する場合、通常、振幅値の誤差は±０．０５％以内、角度誤差は±０．０５°以内になっている。

（２）実例１と実例２とを比較すれば分かるように、ＳＣＡＤＡまたは状態推定断面間の時間間隔を短くすることにより、ダイナミクスのリアルタイム推定の精度を向上することができる。また、実例２から分かるように、通常のＳＣＡＤＡの最大のサンプリング間隔である５秒を用いた場合においても、本発明の方法は、非常に高い推定精度を実現することができる。

（３）実例３から分かるように、ＰＭＵ有りの測定点と推定しようとする点との間の距離が比較的長い場合、誤差はすこし大きくなっているが、その変化があまり著しくない。これは、推定しようとする量の近くにＰＭＵ測定点が存在しない場合、推定しようとする量から離れた同一電気島に存在するその他のノードのＰＭＵ測定を利用しても、ダイナミクスのリアルタイム推定を行うことができることを意味している。

（４）実例４から分かるように、異なる電圧レベルの測定量の間にも、ダイナミクスのリアルタイム推定を行うことができ、しかも、高い推定精度を実現することができる。

（５）実例５から分かるように、位相角のダイナミクスのリアルタイム推定の誤差は約±０．０５°であり、非常に高い精度が実現されている。その理由は、異なるノードの位相角変化量の間に、強い線形関係が存在しているからである。

（６）実例６および実例７から分かるように、電圧量に比べて、有効電力の方は、ダイナミクスのリアルタイム推定の誤差が大きい。その理由は、有効電力の変化量の間に、線形関係が比較的弱く、即ち、厳密な線形関係が存在しないからである。それにしても、上記誤差は、工程応用上の許容可能な誤差範囲内に収まっている。例えば、厳しいダイナミクス推定の状況において、即ち、実例７における低周波振動期間中の非協同グループ間の有効電力量に対するダイナミクス推定の場合においても、誤差を±５％内に維持することができる。

（７）さらに、各実例から分かるように、再帰最小二乗法によるダイナミクスのリアルタイム推定の初期段階において、線形結合係数の初期推定値の誤差が大きいため、ダイナミクスのリアルタイム推定の誤差も相対的に大きい。しかしながら、ダイナミクス推定の進行に伴い、線形結合係数の初期推定値の誤差による影響が小さくなって、ダイナミクスのリアルタイム推定の誤差も著しく改善されている。

Claims (10)

1. ネットワークのパラメータおよび状態マトリックスに依存せず、フェーザ計測ユニットであるＰＭＵを備えた測定点の実測動態データ、および電力系統データ収集と監視制御システムであるＳＣＡＤＡの測定データまたは状態推定データのみに基づいて、ＰＭＵを備えていない低可観測性測定点のダイナミクスをリアルタイムで推定する方法であって、当該ダイナミクスのリアルタイム推定方法は、再帰最小二乗法を利用して、一定期間における準定常状態のデータシーケンス、即ち、連続した複数の断面のＳＣＡＤＡデータまたは状態推定データに基づいて、ＰＭＵを備えていないノードの推定しようとする量の変化量と、ＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合関係を見つけ、そして、当該線形結合関係を利用して、ＰＭＵを備えていないノードの推定しようとする量のダイナミクスをリアルタイムで推定し、当該ダイナミクスのリアルタイム推定方法は、
   ダイナミクスの推定を行う必要がある、ＰＭＵを備えていないノードiを選択し、当該ノードの電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数および発電機回転速度を含む、当該ノードにおけるダイナミクスを推定する必要がある物理量X_iを選択する、第１のステップと、
   各ＰＭＵを備えていないノードi毎に、そのダイナミクスを推定するために使われる、ＰＭＵを備えたノードの集合S_PMU-iを選択し、当該ノードの集合が、ＰＭＵを備えていないノードとの電気的距離が比較的近い複数の、ＰＭＵを備えたノードからなり、ＰＭＵを備えたノードと推定しようとするノードとが同じ電気島に属していることが保証される前提で、前記集合S_PMU-i中のノードを、オフラインで形成された典型的な運転方式におけるノードインピーダンスマトリックスに基づいて選択し、または手動で指定することが可能であり、ＰＭＵを備えたノードとＰＭＵを備えていないノードとの間の最小線路段数に基づいて、オンライン自動探索方式で決定することも可能であり、集合S_PMU-i中のノードの、推定しようとする測定と同じ種類の測定を、ダイナミクス推定用に選択することが可能である、第２のステップと、
   各ＰＭＵを備えていないノードiにおける推定しようとする測定量X_iの変化量ΔX_iと該当のＰＭＵを備えたノードにおける対応の測定量X_jの変化量ΔX_jとの間の線形結合係数b_ij、および再帰最小二乗法パラメータに対して、初期化を行い、但し、j∈S_PMU-iとする、第３のステップと、
   各推定しようとする測定X_i毎に、現在最新のＳＣＡＤＡ測定時刻T_kにおける計測値または最新の状態推定時刻T_kにおける計算値に基づき、さらに、ヒストリー時刻におけるＳＣＡＤＡ計測値または状態推定値を参照して、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔX_i(T_k)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔX_j(T_k)との間の線形結合係数b_ij(T_k)を算出し、但し、j∈S_PMU-iとし、T_kの下付き「k」が第k回のＳＣＡＤＡ測定または第k回の状態推定を示し、循環変数kが0から始まる整数を示す、第４のステップと、
   第４のステップにおいて、T_k時刻におけるＳＣＡＤＡ計測値または状態推定値に基づいて算出された、現在推定しようとする量の変化量ΔX_iとＰＭＵを備えたノードにおける対応の測定量の変化量ΔX_jとの間の線形結合係数b_ij(T_k)に基づいて、以下の式
   

   

   を用いて、次回のＳＣＡＤＡ測定時刻T_k+1または状態推定時刻T_k+1までの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵを備えていないノードの推定しようとする量X_i(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻T_kにおける該当の量X_i(T_k)に対する変化量の推定値ΔX_i(t)を算出し、但し、T_k<t<T_k+1とし、さらに、方程式X_i(t)=ΔX_i(t)+X_i(T_k) に基づいて、当該ＰＭＵを備えていないノードの推定しようとする量の各ＰＭＵ測定時刻における推定値X_i(t)を算出し、但し、T_k<t<T_k+1とする、第５のステップと、
   新しいＳＣＡＤＡ測定データまたは状態推定結果が得られると、第４のステップに戻り、ＰＭＵを備えていないノードの推定しようとする量の変化量とＰＭＵを備えたノードに
   おける対応の測定量の変化量との間の線形結合係数を更新する、第６のステップとを含む、ことを特徴とするダイナミクスのリアルタイム推定方法。
2. 第２のステップにおいて、前記ＰＭＵを備えていないノードiのダイナミクスを推定するために使われる、ＰＭＵを備えたノードの集合S_PMU-iは、選択された４〜６つのノードから構成されている、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
3. 第２のステップにおいて、前記オンライン自動探索方法は、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、ノードiと１段の線路を介して接続、即ち直接接続しているＰＭＵノードを集合S_PMU-iの中に優先に入れてから、前記集合中のノード数が予め決定された数になるまで、ノードiとそれぞれ２段の線路、３段の線路、４段の線路を介して接続しているＰＭＵを備えたノードを、介した線路の段数が小さい方から順に集合S_PMU-iに加える方法である、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
4. 第２のステップにおいて、推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵを備えたノードの電圧振幅値または電流振幅値を、ＰＭＵを備えていないノードの電圧レベルにおける電圧振幅値または電流振幅値に変換する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
5. 第２のステップにおいて、推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するノードの電圧フェーザを基準とした相対位相角に変換し、即ち、相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
6. 第２のステップにおいて、推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、さらに、変圧器の結線方式によって生じた位相シフトに基づいて、ＰＭＵを備えたノードの電圧位相角または電流位相角を、ＰＭＵを備えていないノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
7. 第３のステップにおける初期化の内容は、下記の処理（ａ）、処理（ｂ）、処理（ｃ）、処理（ｄ）をさらに含み、
   処理（ａ）：まず、推定しようとする測定量の変化量ΔX_iの線形結合係数b_ijを求めるという問題の再帰最小二乗法の数学モデルを作成し、再帰最小二乗法を用いて問題を解く場合の一般的な形式は、
   

   

   であり、dが観測された出力ベクトル、wが求めようとするパラメータベクトル、Aが複数のベクトルu_jからなる状態マトリックスであり、また、u_jがパラメータベクトルwと観測量d( j )との関係を表し、vがノイズベクトルであり、推定しようとする測定量の変化量ΔX_iの線形結合係数b_ijを求めるという問題について、前記式において、計算中に必要な変数の定義は、
   

   

   となっており、ΔX_iは、推定しようとする測定X_iの現在のＳＣＡＤＡ測定値と前回のＳＣＡＤＡ測定値または状態推定値との差値、または推定しようとする測定X_iの現在の状態推
   定値と前回の状態推定値との差値であり、各推定しようとする測定をそれぞれ推定するため、実際には、出力ベクトルdは１つの要素しか存在せず、ΔX₁ ΔX₂ … ΔX_j … ΔX_Niは、それぞれ、S_PMU-i中の各ノードの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値であり、b_i1 b_i2 … b_ij … b_iNi は、それぞれ、S_PMU-i中の各ノードの該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔX_iとの間の線形結合係数であり、
   処理（ｂ）：前記推定しようとする測定量の変化量ΔX_iの線形結合係数 b_ijの初期値を1/N_iとすると、再帰最小二乗法で求めようとするパラメータベクトルwの初期値はw_-1=[1/N_i 1/N_i … 1/N_i ]^Tとなり、但し、N_iはＰＭＵを備えていないノードiのS_PMU-i集合中の要素の数であり、下付き「-1」は初期値を表すものであり、
   処理（ｃ）：前記各線形結合係数と対応する、ＰＭＵを備えた測定の変化量に対して、重み付け係数を設定することにより、N_i×N_i次の対角マトリックスとなる重みマトリックスΠが得られ、その対角要素は各重み付け係数からなり、ＰＭＵを備えた測定の変化量の測定品質が高いほど、また、そのノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きく、再帰最小二乗法の計算プロセスにおける中間プロセスマトリックス変数Pの初期値P_-1をΠとし、ＰＭＵを備えた測定の変化量の測定品質が不明な場合、各重み付け係数を同じ値の０．１とし、この場合、
   

   

   となり、但し、下付「-1」は初期値を表すものであり、
   処理（ｄ）： 望ましくは、初期化忘却因子λを０．９９とする、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
8. 第４のステップにおいて、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔX_i(T_k)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔX_j(T_k)との間の線形結合係数b_ij(T_k)を求める計算には、下記のステップ（ａ）、ステップ（ｂ）、ステップ（ｂ） 、ステップ（ｃ）を含み、
   ステップ（ａ）：ノードiの推定しようとする測定の現在の第k時間ステップ、即ちT_k時刻におけるＳＣＡＤＡ測定値または状態推定値の、前回のサンプリング時刻におけるＳＣＡＤＡ測定値または状態推定値である第k-1時間ステップにおける対応値に対する変化量を算出することにより、第k時間ステップの出力ベクトルd_kが得られ、即ち、
   

   

   

   ステップ（ｃ）：次の３つの回帰公式に基づいて、第k時間ステップにおける再帰最小二乘推定法の中間プロセスマトリックスg_k、パラメータベクトルマトリックスw_k=[ b_i1(T_k) b_i1(T_k) … b_iNi(T_k) ]^T、中間プロセスマトリックスP_kを順に算出し、循環変数kは第0時間ステップから始まり、第0時間ステップの計算に使われたP_-1およびw_-1は、第３のステップにおける初期化処理プロセスによって得られるものであり、第k時間ステップで得られたw_kによって、第k時間ステップの各線
   

   

   形組み合わせ係数b_ij(T_k)が得られる、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
9. 第５のステップにおいて、次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間に、即時に更新されたＳＣＡＤＡスイッチング量が、電力ネットワークに短絡事故または断線事故が発生したことを示した場合、さらに、ＰＭＵを備えていないノードiまたはＰＭＵを備えたノードj（j∈S_PMU-i）に直接接続されている線路に、短路または断線が発生した場合、下記の処理（ａ）、処理（ｂ）および処理（ｃ）を実施し、
   処理（ａ）：ＰＭＵを備えたノードjが系統から孤立された場合、それと対応した線形結合係数b_ijをその他のＰＭＵを備えたノードに平均配分し、孤立されたＰＭＵを備えたノードjの線形結合係数b_ijを0とし、
   処理（ｂ）：ＰＭＵを備えていないノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において当該ノードのダイナミクスを推定することが不可能であることを示し、孤立されたＰＭＵを備えていないノードiが再び系統に接続されたら、孤立される前の線形結合係数を初期係数b_ijとして、該当の測定のダイナミクスを推定し、
   処理（ｃ）：その他の状況、即ち、短絡または断線によってＰＭＵを備えていないノードiまたはＰＭＵを備えたノードjが系統から孤立されたという状況にならなかった場合、特殊な処理を行わず、次回の、即ち第k+1時間ステップのＳＣＡＤＡアナログ量測定データまたは状態推定データの更新によって、線形結合係数b_ijを補正して、ネットワーク構成の変化を表す、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
10. 第５のステップにおいて、推定された相対位相角値を-180°〜+180°の範囲内に換算し、または、推定された相対位相角値を絶対位相角値に変換し、さらに、-180°〜+180°の範囲内に換算する、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。

Images