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

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

Info

Publication number
JP5486686B2
JP5486686B2 JP2012530072A JP2012530072A JP5486686B2 JP 5486686 B2 JP5486686 B2 JP 5486686B2 JP 2012530072 A JP2012530072 A JP 2012530072A JP 2012530072 A JP2012530072 A JP 2012530072A JP 5486686 B2 JP5486686 B2 JP 5486686B2
Authority
JP
Japan
Prior art keywords
pmu
node
measurement
estimated
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2012530072A
Other languages
English (en)
Other versions
JP2013506389A (ja
Inventor
ドアン・ガン
ウー・ジンタウ
スウン・シヤオイエン
ヤン・ドン
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Sifang Automation Co Ltd
Beijing Sifang Engineering Co Ltd
Original Assignee
Beijing Sifang Automation Co Ltd
Beijing Sifang Engineering Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Sifang Automation Co Ltd, Beijing Sifang Engineering Co Ltd filed Critical Beijing Sifang Automation Co Ltd
Publication of JP2013506389A publication Critical patent/JP2013506389A/ja
Application granted granted Critical
Publication of JP5486686B2 publication Critical patent/JP5486686B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J13/00—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network
    • H02J13/00002—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network characterised by monitoring
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J13/00—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network
    • H02J13/00001—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network characterised by the display of information or by user interaction, e.g. supervisory control and data acquisition systems [SCADA] or graphical user interfaces [GUI]
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J13/00—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network
    • H02J13/00032—Systems characterised by the controlled or operated power network elements or equipment, the power network elements or equipment not otherwise provided for
    • H02J13/00034—Systems characterised by the controlled or operated power network elements or equipment, the power network elements or equipment not otherwise provided for the elements or equipment being or involving an electric power substation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R19/00—Arrangements for measuring currents or voltages or for indicating presence or sign thereof
    • G01R19/25—Arrangements for measuring currents or voltages or for indicating presence or sign thereof using digital measurement techniques
    • G01R19/2513—Arrangements for monitoring electric power systems, e.g. power lines or loads; Logging
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
    • Y02E40/70—Smart grids as climate change mitigation technology in the energy generation sector
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/22—Flexible AC transmission systems [FACTS] or power factor or reactive power compensating or correcting units
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/30—State monitoring, e.g. fault, temperature monitoring, insulator monitoring, corona discharge
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/40—Display of information, e.g. of data or controls

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Human Computer Interaction (AREA)
  • Supply And Distribution Of Alternating Current (AREA)
  • Remote Monitoring And Control Of Power-Distribution Networks (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、電力系統の広域測定技術分野に属し、特に、ＰＭＵを備えていない測定点の測定量のダイナミクスを推定する方法に関する。
従来、電力系統コントロールセンターの情報は、主に、ＳＣＡＤＡシステム（データ収集および監視制御システム）から収集されたものであり、電圧、電流、電力、周波数などのリモート測定のサンプリング間隔が１〜５秒であり、しかも、タイムスケールがないため、電力ネットワークの準定常状態の情報しか反映することができない。近年来、グローバル・ポジショニング・システムであるＧＰＳに基づく同期フェーザ測定技術の電力系統への応用が益々広がり、当該技術に基づくコアデバイスは、各変電所および発電所に配置されるフェーザ測定ユニット、即ちＰＭＵである。ＰＭＵは、各変電所・発電所における、同期タイムスケールがマークされた電圧フェーザ、電流フェーザ、電力、周波数、発電機回転速度などの情報を、１０ｍｓ〜３０ｍｓのサンプリング間隔でコントロールセンターのマスターステーションに送信する。これによって、コントロールセンターにおいて、電力系統の動的挙動に対するオンライン集中監視が実現される。しかしながら、現在のＰＭＵ装置の価額がかなり高く、尚且つ、ＰＭＵのために大量な通信資源を装備する必要があるため、現段階および今後相当長期間にわたって、一部の変電所・発電所にしかＰＭＵを取り付けることができない。そのため、現在、電力ネットワークのコントロールセンターでは、ＰＭＵに基づく広域測定システム（ＷＡＭＳ）において、電力ネットワークの動的挙動に対する可観測性のすべての要求を満たすことがまだ実現されていない。
ＰＭＵを取り付けておらず（または電力会社に管轄されておらず）、しかも完全可観測性を有していない変電所・発電所またはノードのダイナミクスに対するリアルタイム監視を実現するために、補助オペレーターがオンラインで安全・安定対策を実施する。そのため、既存のオンライン測定情報を利用して、ＰＭＵを備えていないノード（以下、ＰＭＵ無しのノードと称することもある）のダイナミクスをリアルタイムで推定することは、非常に重要かつ実用的なことである。現在、様々な文献に開示された方法の中に、電力ネットワークにおける、ＰＭＵを備えたノード（以下、ＰＭＵ有りのノードと称することもある）の電圧フェーザとＰＭＵ無しのノードの電圧フェーザとの間の感度関係、および該当のＳＣＡＤＡ測定情報に基づいて、ＰＭＵ有りのノードの電圧フェーザの振幅値と位相角のダイナミクスによって、ＰＭＵ無しのノードの電圧フェーザの振幅値と位相角のダイナミクスを推定する方法が存在する。但し、ＰＭＵ有りのノードの電圧フェーザ（振幅値と位相角）とＰＭＵ無しのノードの電圧フェーザ（振幅値と位相角）との間の感度関係は、状態マトリックス、即ち潮流方程式のヤコビ行列によって導出される。このような方法では、以下の３つの問題点が存在する。
（１）原始ヤコビ行列の各要素が、ノード電圧の振幅値、位相角およびノードアドミタンスと関連しているため、潮流およびネットワーク状態が変化すると、それに伴ってヤコビ行列を常に更新する必要がある。しかしながら、ヤコビ行列更新の計算量が膨大であり、オンラインで１０ｍｓレベルのダイナミクスのデータを生成することは適切ではない。それに対して、電圧振幅値および位相角の変化を無視し、さらに、電力ネットワーク構成要素の抵抗を無視するなど一連の簡略化措置を取ることにより、ヤコビ行列の要素をノードのアドミタンス行列の虚部の要素に取り替えて、ＰＭＵ無しのノードとＰＭＵ有りのノードとの間の近似線形相関係数を求める方法が、ある文献に開示されている。なお、これらの相関係数は、ネットワークの構造が変わった時のみ更新され、その他の状況における潮流変化は、ノード間の相関係数の更新を引き起こさない。これによって、オンライン計算
の作業量は大幅に減少されるが、ダイナミクスのオンライン推定の精度も一定程度下がる。
（２）上記したように、簡略化されたヤコビ行列をアドミタンス行列から得ることはできるが、アドミタンス行列は、線路、変圧器、発電機などのネットワーク構成要素のアドミタンスパラメータに依存している。通常、これらのアドミタンスパラメータは、メーカが提供した定格値または典型的な値から得られたものである。これらの数値と実際のネットワーク構成要素のパラメータとの間は、時々、大きな誤差が存在する。しかも、これらの構成要素のパラメータは、温度、周波数などの運転環境または運転状態の変化に伴って変化することもある。これらのパラメータが精確でなければ、ＰＭＵ有りのノードとＰＭＵ無しのノードとの測定量の変化量間の線形相関関係を精確に計算することができず、さらに、ＰＭＵ無しのノードの測定量のダイナミクスを精確に推定することができない。また、電力ネットワークの構成要素が非常に多いので、上記パラメータが精確か否かについて、一つずつ検証してオンライン補正を行うことが不可能なため、ダイナミクスの推定が不精確になっても補正されない。
（３）既存の文献に開示されたダイナミクスの推定は、電圧フェーザのみに限られており、その他の測定、特に電力に対する処理方法に言及していない。
また、既存の文献において、ＰＭＵおよびＳＣＡＤＡデータに基づいて混合状態推定または混合動態状態推定を行うといった内容は多いが、いずれもＰＭＵデータを利用して、現在または将来のＳＣＡＤＡ測定断面の状態を推定する精度を向上するものであり、ＰＭＵ無しのノードの測定に対するダイナミクスの推定を行うという内容に関連していないと言える。
したがって、上記説明から分かるように、現在、ＰＭＵ無しの測定点の測定に対するダイナミクスの推定を行うための、精確で実用的な方法はまだ存在していない。このような状況において、本発明は、ネットワークのパラメータおよび状態マトリックスの更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、フェーザ測定ユニットであるＰＭＵを備えた測定点の実測動態データおよび電力系統データ収集と監視制御システムであるＳＣＡＤＡの測定データまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のダイナミクスをリアルタイムで推定する方法を提供する。本発明の方法は、線形化理論を理論基礎として、現在の断面とヒストリー断面のＳＣＡＤＡ（または状態推定）データに基づいて、最小二乗法を利用して、ＰＭＵ有りのノードの測定量の変化量とＰＭＵ無しのノードの対応測定量の変化量との間の線形結合関係を見つけ出し、そして、ＰＭＵ測定時間断面において、ＰＭＵ有りのノードの測定情報を用いて、ＰＭＵ無しのノードの該当する測定情報を推定する。ダイナミクス推定の全体のプロセスは、完全に実測データに依拠しており、電力ネットワークの構成要素のパラメータおよびヤコビ行列または状態マトリックスの更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がない。測定量の変化量間の線形関係が比較的強い場合、当該方法によるダイナミクス推定の誤差は非常に小さい。例えば、電圧フェーザのダイナミクスに対する推定の場合、通常、振幅値の誤差は±０．０５％以内、角度誤差は±０．０５°以内になっている。一方、測定量の変化量間の線形関係が比較的弱い場合、ダイナミクス推定の誤差は比較的大きいが、工程応用上の許容誤差範囲内に止まる。例えば、通常、低周波振動期間における電力ネットワークの有効電力に対するダイナミクス推定の誤差は、約±５％内に収まることが可能である。
ＰＭＵ無しの測定点のダイナミクスを推定する方法の精確性が、ネットワークのパラメータおよびヤコビ行列のリアルタイム更新に依存し、しかも、推定が電圧フェーザのみに限られているという従来技術の短所を克服するために、本発明は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、ＰＭＵ有りの測定点
の実測動態データおよびＳＣＡＤＡ測定データまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のダイナミクスをリアルタイムで推定する方法を開示する。当該方法では、ＰＭＵ無しの測定点は、完全可観測性を有する必要がない。本発明は、以下の具体的な技術案を用いる。
本発明に係るダイナミクスのリアルタイム推定方法は、一定期間における準定常状態のデータシーケンス、即ち、連続した複数の断面のＳＣＡＤＡデータまたは状態推定データに基づいて、ＰＭＵを備えていないノードの推定しようとする量の変化量と、ＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合関係を見つけ、そして、当該線形結合関係を利用して、ＰＭＵを備えていないノードの推定しようとする量のダイナミクスをリアルタイムで推定する。本発明を実現するための典型的なフローチャートは図１Ａと図１Ｂに示すようになっており、具体的な実施ステップは、以下の通りである。
第１のステップ：ダイナミクスの推定を行う必要があり、ＰＭＵを備えていないノードiを選択し、当該ノードにおける該当の構成要素の電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数および発電機回転速度などの、当該ノードにおけるダイナミクスを推定する必要がある物理量Xiを選択する。
第２のステップ：各ＰＭＵを備えていないノードi毎に、そのダイナミクスを推定するために使われ、ＰＭＵを備えたノードの集合SPMU-iを選択し、当該ノードの集合は、通常、ＰＭＵを備えていないノードとの電気的距離が比較的近い複数の、ＰＭＵを備えたノードからなり、実用上、４〜６つのノードとすることが好ましい。本方法において、電気的距離に対する要求は厳しくなく、ＰＭＵを備えたノードと推定しようとするノードとが同じ電気島（electrical island）に属していることが保証される前提で、集合SPMU-i中のノードを、オフラインで形成された典型的な運転方式におけるノードインピーダンスマトリックスに基づいて選択し、または手動で指定することができ、ＰＭＵを備えたノードとＰＭＵを備えていないノードとの間の最小線路段数に基づいて、オンライン自動探索方式で決定してもよい。具体的なステップとして、まず、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、ノードiと１段の線路を介して接続、即ち直接接続しているＰＭＵノードを集合SPMU-iの中に優先に入れてから、集合中のノード数が予め決定された数になるまで、ノードiとそれぞれ２段の線路、３段の線路、４段の線路を介して接続しているＰＭＵを備えたノードを、介した線路の段数が小さい方から順に集合SPMU-iに加える。集合SPMU-i中のノードの、推定しようとする測定と同じ種類の測定をダイナミクス推定用に選択することができる。
第２のステップにおいて、次のデータ予備処理プロセスがさらに含まれている。
（ａ）推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵを備えたノードの電圧振幅値または電流振幅値を、ＰＭＵを備えていないノードの電圧レベルにおける電圧振幅値または電流振幅値に変換しなければならない。
（ｂ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するあるノードの電圧フェーザを基準とした相対位相角に変換し、当該相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算するようにしなければならない。
（ｃ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、変圧器の
結線方式によって生じた位相シフトに基づいて、ＰＭＵを備えたノードの電圧位相角または電流位相角を、ＰＭＵを備えていないノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算するようにしなければならない。
第３のステップ：各ＰＭＵを備えていないノードiにおける推定しようとする測定量Xiの変化量ΔXiと該当のＰＭＵを備えたノードにおける対応の測定量Xjの変化量ΔXjとの間の線形結合係数bij、および再帰最小二乗法パラメータに対して、初期化を行い、但し、j∈SPMU-iとする。初期化の具体的な方法は以下の通りである。
（ａ）まず、推定しようとする測定量の変化量ΔXiの線形結合係数bijを求めるという問題の再帰最小二乗法の数学モデルを作成する。再帰最小二乗法を用いて問題を解く場合の一般的な形式は、
Figure 0005486686
であり、但し、dは観測された出力ベクトル、wは求めようとするパラメータベクトル、Aは複数のベクトルujからなる状態マトリックスであり、ujはパラメータベクトルwと観測量d( j )との関係を表すものであり、vはノイズベクトルである。推定しようとする測定量の変化量ΔXiの線形結合係数bijを求めるという問題について、上式における計算中に必要な変数の定義は、
Figure 0005486686
となっており、ΔXiは、推定しようとする測定量Xiの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との間の差値であり、各推定しようとする測定をそれぞれ推定するため、実際には、出力ベクトルdは１つの要素しかない。ΔX1 ΔX2 …ΔXj …ΔXN は、それぞれ、SPMU-i中の各ノードの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値であり、bi1 bi2 … bij … biNi は、それぞれ、SPMU-i中の各ノードの該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔXiとの間の線形結合係数である。
（ｂ）前記推定しようとする測定量の変化量ΔXiの線形結合係数 bijの初期値を、通常1/Niとすると、再帰最小二乗法で求めようとするパラメータベクトルwの初期値はw-1=[1/Ni 1/Ni … 1/Ni ]Tとなり、但し、NiはＰＭＵを備えていないノードiのSPMU-i集合中の
要素の数であり、下付き「-1」は初期値を表す。
（ｃ）前記各線形結合係数と対応する、ＰＭＵを備えた測定の変化量に対して、重み付け係数を設定することにより、Ni×Ni次の対角マトリックスである重みマトリックスΠが得られ、その対角要素は、各重み付け係数からなる。ＰＭＵを備えた測定の変化量の測定品質が高いほど、また、そのノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きい。再帰最小二乗法の計算プロセスにおける中間プロセスマトリックス変数Pの初期値P-1をΠとし、ＰＭＵを備えた測定の変化量の測定品質が不明な場合、各重み付け係数を同じ値の0.1とし、この場合、
Figure 0005486686
となり、但し、下付き「-1」は初期値を表す。
（ｄ）望ましくは、初期化忘却因子λを０．９９とする。
第４のステップ：各推定しようとする測定量Xi毎に、現在最新のＳＣＡＤＡ測定時刻Tkにおける計測値または最新の状態推定時刻Tkにおける計算値、およびヒストリー時刻におけるＳＣＡＤＡ計測値または状態推定値に基づき、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔXi(Tk)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔXj(Tk)との間の線形結合係数bij(Tk)を算出し、但し、j∈SPMU-iとし、Tkの下付き「k」は第k回のＳＣＡＤＡ測定または第k回の状態推定を示し、循環変数kは0から始まる整数を示す。具体的には、以下のステップ（ａ）〜（ｃ）を含む。
（ａ）ノードiの推定しようとする測定の現在の第k時間ステップ、即ちTk時刻におけるＳＣＡＤＡ測定値または状態推定値の、前回のサンプリング時刻におけるＳＣＡＤＡ測定値または状態推定値である第k-1時間ステップにおける対応値に対する変化量を算出することにより、下記第k時間ステップの出力ベクトルdkが得られる。
Figure 0005486686
（ｂ）次の式により、集合SPMU-i中の各ノードの、推定しようとする測定と対応する各測定の現在の第k時間ステップにおけるＳＣＡＤＡ測定値または状態推定値の、第k-1時間ステップにおける対応値に対する変化量のベクトルukを算出する。
Figure 0005486686
（ｃ）次の３つの回帰公式に基づいて、第k時間ステップにおける再帰最小二乘推定法
の中間プロセスマトリックスgk、パラメータベクトルマトリックスwk=[ bi1(Tk) bi1(Tk) … biNi(Tk) ]T、中間プロセスマトリックスPkを順に算出する。循環変数kは第0時間ステップから始まり、第0時間ステップの計算に使われたP-1およびw-1は、第３のステップにおける初期化処理プロセスによって得られるものである。第k時間ステップにおける各線形組み合わせ係数bij(Tk)は、第k時間ステップで得られたwkによって得られる。
Figure 0005486686
第５のステップ：第４のステップにおいて、Tk時刻におけるＳＣＡＤＡ計測値または状態推定値に基づいて算出された、現在推定しようとする量の変化量ΔXiとＰＭＵを備えたノードにおける対応の測定量の変化量ΔXjとの間の線形結合係数bij(Tk)に基づいて、以下の式
Figure 0005486686
を用いて、次回のＳＣＡＤＡ測定時刻Tk+1または状態推定時刻Tk+1までの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵを備えていないノードの推定しようとする量Xi(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻Tkにおける該当の量Xi(Tk)に対する変化量の推定値ΔXi(t)を算出し、但し、Tk<t<Tk+1とする。
さらに、方程式Xi(t)=ΔXi(t)+Xi(Tk) に基づいて、当該ＰＭＵを備えていないノードの推定しようとする量の各ＰＭＵ測定時刻における推定値Xi(t)を算出し、但し、Tk<t<Tk+1とする。なお、推定された相対位相角の値を-180°〜+180°の範囲内に換算し、または、推定された相対位相角の値を絶対位相角の値に変換してから、-180°〜+180°の範囲内に換算するようにしなければならない。
ＳＣＡＤＡ測定または状態推定値が存在する時刻において、ＰＭＵ無しのノートの推定しようとする量は、ＳＣＡＤＡ測定値または状態推定値とする。
さらに、第５のステップにおいて、次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間に、電力ネットワークに短絡または断線事故が発生する可能性があり、事故が発生した場合、ＳＣＡＤＡの該当スイッチング量が即時に更新されるが、アナログ量は、通常、即時に更新されない。この場合、ネットワーク構成の変化により、ヤコビ行列またはノードアドミタンス行列が変わるため、ＰＭＵを備えた各ノードによる、ＰＭＵを備えていないノードの推定量の変化への影響は大きく変化する。しかしながら、ＳＣＡＤＡのアナログ測定データまたは状態推定データの更新が遅いため、新しい係数を直ぐに生成することが難しい。このような推定上の誤差を減らすために、トポロジー解析によって、ＰＭＵを備えていないノードiまたはＰＭＵを備えたノードj（j∈SPMU-i）に直接接続されている線路に、短路または断線が発生した場合、以下の処理（ａ）〜（ｃ）を実施する。即ち、
（ａ）ＰＭＵを備えたノードjが系統から孤立された場合、それと対応した線形結合係
数bijをその他のＰＭＵを備えたノードに平均配分し、孤立されたＰＭＵを備えたノードjの線形結合係数bijを0とし、
（ｂ）ＰＭＵを備えていないノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において当該ノードのダイナミクスを推定することが不可能であることを示し、孤立されたＰＭＵを備えていないノードiが再び系統に接続されたら、孤立される前の線形結合係数を初期係数bijとして、該当の測定のダイナミクスを推定し、
（ｃ）その他の状況、即ち、短絡または断線によってＰＭＵを備えていないノードiまたはＰＭＵを備えたノードjが系統から孤立されたという状況にならなかった場合、特殊な処理を行わず、次回の、即ち第k+1時間ステップのＳＣＡＤＡアナログ測定データまたは状態推定データの更新によって、線形結合係数bijを補正して、ネットワーク構成の変化を表す。
第６のステップ：新しいＳＣＡＤＡ測定データまたは状態推定結果が得られると、第４のステップに戻り、ＰＭＵを備えていないノードの推定しようとする量の変化量とＰＭＵを備えたノードにおける対応の測定量の変化量との間の線形結合係数を更新する。
本発明は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、ＰＭＵ有りの測定点の実測動態データおよびＳＣＡＤＡデータまたは状態推定データのみに基づいて、ＰＭＵ無しの測定点のＰＭＵ測定期間におけるダイナミクスのリアルタイム推定を実現することができ、従来のダイナミクス推定方法における、精確性が電力ネットワークパラメータおよびヤコビ行列の精確性に依存しているため、計算量が大きく、リアルタイム性が低く、精確性が確保できないといった問題点を克服することができる。本発明の方法を利用する場合、測定量の変化量間の線形関係が比較的強い時に、ダイナミクス推定の誤差は非常に小さい。例えば、電圧フェーザに対するダイナミクス推定の場合、通常、その振幅値の誤差は±０．０５％以内、角度の誤差は±０．０５°以内である。一方、測定量の変化量間の線形関係が比較的弱い時に、ダイナミクス推定の誤差は相対的に大きいが、工程応用上の許容可能な誤差範囲内に止まる。例えば、低周波振動期間における電力ネットワークの有効電力に対するダイナミクスの推定の場合、通常、その誤差を約±５％以内に保つことができる。なお、ＳＣＡＤＡまたは状態推定データ断面間の時間間隔の短縮に伴って、上述したダイナミクス推定の誤差がさらに小さくなる。
状態マトリックスに依存せずにＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法のフローチャートである。 状態マトリックスに依存せずにＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法のフローチャートである。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（近くのＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が１秒である）。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（近くのＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が５秒である）。 光谷ステーション夙光線の電圧振幅値ダイナミクスの推定値と実測値との相対誤差グラフである（離れたＰＭＵステーションを利用して推定され、ＳＣＡＤＡ周期が１秒である）。 ５００ＫＶ側測定によって光谷ステーション４号変圧器２２０ＫＶ母線の電圧振幅値ダイナミクスを推定した場合の相対誤差グラフである。 光谷ステーション夙光線の電圧位相角ダイナミクスの推定値と実測値との絶対誤差グラフである。 光谷ステーション夙光線の有効電力ダイナミクスの推定値と実測値との相対誤差グラフである。 ある低周波振動の前および振動中に、非協同領域の振動電力を用いて万龍１線の有効電力ダイナミクスを推定した場合の相対誤差グラフである。
本発明に開示される方法は、ネットワークパラメータおよび状態マトリックス（またはヤコビ行列）の更新に依存せず、しかも、ＰＭＵ無しの測定点が完全可観測性を有する必要がなく、ＰＭＵ有りの測定点の実測動態データおよびＳＣＡＤＡ測定データまたは状態推定データのみに基づいて、ＰＭＵ無しのノードの測定量のダイナミクスに対するリアルタイム推定を実現することができる。本発明は、ある測定量の変化量を、当該測定量と相関する複数の測定量の変化量の線形結合から得ることができるという発想に基づくものである。即ち、
Figure 0005486686
上式において、Xiはノードiにおける、ダイナミクスが推定される１つの測定量であり、電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数または発電機回転速度などであってもよい。ΔXi(t)は、時刻tにおけるXiの変化量を表すものである。集合SPMU-iは、例えば、変電所および発電所のトポロジーノードのような、ノードiとの電気的距離が近いＰＭＵ有りのノードからなる。SPMU-i中のノードに対して、ノードiとの電気的距離が近いことを要求しているが、この要求は厳格な要求ではない。後述の計算例から分かるように、ノードiとの電気的距離が比較的遠いが、ノードiと同一電気島内に存在するノードを利用しても、同様に、ノードiのダイナミクスに対する推定を実現することができる。この場合、誤差は大きくなるが、通常の工程応用上の誤差許容範囲内に収まる。Xjは集合SPMU-i中のノードjにおける、Xiと同じ性質の測定量であり、bijは相関する測定量の変化量の線形結合係数である。
本方法において、集合SPMU-iの大きさに関する厳密な規定はない。式（１）から分かるように、SPMU-i中の要素が1になってはいけない。さもなければ、推定されたノードのダイナミクスの形状は、ＰＭＵ有りのノードのダイナミクスと完全に一致してしまう。様々な要因による、推定されるノードのダイナミクスの形状および大きさへの影響を反映するために、集合SPMU-i中の要素として、数が多く、分布範囲が広いことが望ましい。しかしながら、集合SPMU-i中の要素の増加に伴って、ダイナミクスのリアルタイム推定の計算量も増加するため、ダイナミクス推定のリアルタイム性に影響を与える可能性がある。実際の工程応用において、集合中の要素の数を４〜６とすることが最も適切であると実験によって判明されている。
また、式（１）に示した線形結合関係は、主に経験に基づいた公式であり、厳密な数学関係に依存するものではないため、測定量の変化量間の線形関係が著しい電圧フェーザのダイナミクス推定に応用することができ、電力、周波数など変化量間の線形関係が著しくない測定量のダイナミクス推定に応用することもできる。ユーザは、線形化する具体的な数学モデルに関わる必要がなく、再帰最小二乗法を用いてヒストリーデータ学習を行うことによって、線形化係数を求めることができる。実験で証明されたように、本発明の方法
によれば、ダイナミクス推定に関与する測定量の変化量の間に厳密な線形関係が存在する場合、動態推定結果の精度は非常に高く、例えば、電圧振幅値のダイナミクス推定の誤差は、通常、±０．０５％以内であり、電圧位相角度のダイナミクス推定の誤差は、通常、±０．０５°以内である。一方、ダイナミクス推定に関与する測定量の変化量間の線形関係が比較的に弱い場合、ダイナミクス推定の誤差は比較的に大きいが、工程応用上の誤差許容範囲内に収まり、例えば、低周波振動期間における電力ネットワークの有効電力のダイナミクスを推定する場合でも、通常、誤差を約±５％内に維持することができる。なお、ＳＣＡＤＡまたは状態推定データ断面間の時間間隔の短縮に伴い、例えば、５秒から１秒に変更すると、上記ダイナミクス推定の誤差はさらに小さくなる。
式（１）および上述した原理に基づいて、連続した複数の断面のＳＣＡＤＡ測定データまたは状態推定データを利用して、ＰＭＵ無しのノードの推定しようとする測定量の変化量とＰＭＵ有りのノードにおける対応の測定量の変化量との間の線形結合関係を求める。そして、当該線形結合関係を利用して、ＰＭＵ無しのノードの推定しようとする測定量の各ＰＭＵ測定時間断面における該当値を推定することにより、ＰＭＵ無しのノードの測定量のダイナミクスに対するリアルタイム推定を実現する。
準定常状態のデータシーケンスのソースに関して、電圧振幅値、電流振幅値、有効電力、無効電力および周波数などの電気量のみに対してダイナミクスを推定する場合は、ＳＣＡＤＡデータのみを使用すればよく、通常、そのデータ断面間の間隔は１〜５秒であり、さらに、電圧位相角、電流位相角のダイナミクスを推定する必要がある場合では、状態推定の結果を使用しなければならない。その理由は、現在、ＳＣＡＤＡが位相角の測定を行っていないからである。状態推定データ断面間の時間間隔は、状態推定プログラムの計算効率および電力ネットワークの規模によって異なり、ダイナミクスのリアルタイム推定の応用において、通常、５秒に設定することができる。また、現在、一部のＳＣＡＤＡおよびＰＭＵ／ＷＡＭＳ（広域測定システム）は、発電機の回転速度など非電気量の収集もサポートしている。この場合では、本発明の方法を用いて、ＰＭＵ無しの測定点における非電気量に対するダイナミクス推定を実現することができる。
式（１）に対して、再帰最小二乗法を用いて係数bijを求めるようにするには、次の理由がある。第１、再帰の方法によれば、新しいＳＣＡＤＡまたは状態推定断面を得た直後に当該時間断面の係数bijを求めることができるため、計算のリアルタイム性を維持することができる。第２、式（１）は、以下に示す、再帰最小二乗法で問題を解く場合の一般的な形式に適合している。
Figure 0005486686
式（２）において、dは観測された出力ベクトル、wは求めようとするパラメータベクトル、Aは複数のベクトルujからなる状態マトリックスであり、ujはパラメータベクトルwと観測量d( j )との関係を表すものであり、vはノイズベクトルである。したがって、式（１）と式（２）とは、以下の対応関係を有している。
Figure 0005486686
但し、ΔXiは、推定しようとする測定量Xiの現在のＳＣＡＤＡ測定値（または状態推定値）と前回のＳＣＡＤＡ測定値（または状態推定値）との差値を表すものであり、各推定しようとする測定をそれぞれ推定するため、実際に、出力ベクトルdは1つの要素しかない
。ΔX1 ΔX2 … ΔXj … ΔXNi は、それぞれ、SPMU-i中の各ノードにおける現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値である。
bi1 bi2 … bij … biNi は、それぞれ、SPMU-i中の各ノードにおける該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔXiとの間の線形結合係数である。式（２）は再帰最小二乗法で問題を解く場合の一般的な形式に適合している。そのため、再帰最小二乗法による汎用解法によって、係数bijを便利に算出することができる。
実際のシステムにおいて、本発明の原理に基づくダイナミクスリアルタイム推定ソフトウェアは、以下のような典型的なハードウェア環境で実行することが可能である。コントロールセンターの広域測定システム（ＷＡＭＳ）のハイレベル応用サーバまたはリアルタイムデータサーバにおいて、ダイナミクスリアルタイム推定ソフトウェアを実行することができる。ダイナミクスリアルタイム推定ソフトウェアは、ＳＣＡＤＡ／ＥＭＳシステムのリアルタイムデータベースから、現在の電力ネットワークのＳＣＡＤＡ測定値または状態推定値をオンラインで獲得し（推定しようとする量およびそれと相関するＰＭＵ有りのノードにおける対応の測定量のＳＣＡＤＡ測定値または状態推定値のみを獲得してもよい）、そして、再帰最小二乗法を利用して、現在のＳＣＡＤＡまたは状態推定時間断面に対応する各推定しようとする測定量の各線形結合係数を計算する。次回のＳＣＡＤＡ測定値または状態推定値を受信するまで、ダイナミクスリアルタイム推定ソフトウェアは、現在の線形結合係数およびＷＡＭＳリアルタイムデータサーバから獲得した関連のＰＭＵ有りのノードにおける対応の動態測定値を利用して、ＰＭＵ無しのノードの推定しようとする測定量の各ＰＭＵ測定時間断面における該当値を推定し、そして、各ＰＭＵ測定時間断面における推定結果をＷＡＭＳリアルタイムデータサーバに保存し、さらに、ＰＭＵ無しの実測値であることを示すマークを付ける。
次に、実際のデータを用いて、本発明の技術案の具体的なステップをさらに説明する。図１Ａ、図１Ｂは、本発明の方法をソフトウェアで実現するためのフローチャートを示している。本発明の方法の有効性を検証するために、以下の計算例において、実際にＰＭＵで測定することができる測定量を、ダイナミクスを推定しようとする測定量とし、そして、推定したダイナミクスと実際に測定したダイナミクスとを比較することにより、本発明の方法の有効性を検証するようにしている。ネットワークパラメータおよび状態マトリックスに依存せずに、ＰＭＵ無しの低可観測性測定点のダイナミクスをリアルタイムで推定するアルゴリズムは、以下のステップによって実現される。
（１）ダイナミクスを推定する必要があるＰＭＵ無しのノード、および当該ノードにおける、ダイナミクスを推定する必要がある物理量Xiを選択し、当該物理量Xiは、例えば、ある構成要素の電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効、無効、周波数または発電機回転速度などである。
（２）各ＰＭＵ無しのノードi毎に、そのダイナミクスを推定するために使われる、ＰＭＵ有りのノードの集合SPMU-iを選択する。当該ノードの集合は、通常、ＰＭＵ無しのノードとの電気的距離が比較的近い複数の、ＰＭＵ有りのノードからなり、手動指定の方法または自動探索の方法によって決定される。自動探索方法の１つとして、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、１段の線路を介してノードiと接続、即ち直接接続しているＰＭＵノードを集合SPMU-iに優先に入れて、その後、集合中のノードの数が４〜６（状況によって集合の大きさを選ぶ）になるまで、それぞれ、２段の線路、３段の線路、４段の線路などを介してノードiと接続しているＰＭＵ有りのノードを、介した線路の段数が小さいほうから順に集合SPMU-iに加える。集合SPMU-i中のノードの、推定しようとする測定量と同じ種類の測定をダイナミクスの推定用に選択することができる。推定される量によって、原始測定データに対して、以下の予備処理を行う必要がある。
（２．ａ） 推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵ有りのノードとＰＭＵ無しのノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵ有りのノードの電圧振幅値または電流振幅値を、ＰＭＵ無しのノードの電圧レベルにおける電圧振幅値または電流振幅値に変換しなければならない。
（２．ｂ） 推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するあるノードの電圧フェーザを基準とした相対位相角に変換し、即ち、相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算するようにしなければならない。
（２．ｃ）推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵ有りのノードとＰＭＵ無しのノードとは電圧レベルが異なっているならば、変圧器の結線方式で生じた位相シフトに基づいて、ＰＭＵ有りのノードの電圧位相角または電流位相角を、ＰＭＵ無しのノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算するようにしなければならない。
（３）各ＰＭＵ無しのノードiにおける推定しようとする測定の係数bijおよび再帰最小二乗法のパラメータに対して、初期化を行う。
（３．ａ） 各推定しようとする測定の係数bijの初期値は、1/Niとすることができる。NiはＰＭＵ無しのノードiのSPMU-i集合中の要素の数である。同様に、再帰最小二乗法の未知数ベクトルwに対する初期化を行って、w-1=[1/Ni 1/Ni … 1/Ni ]Tが得られ、下付き「-1」は初期値を表すものである。
（３．ｂ） 各線形結合係数に対応したＰＭＵ有りの測定の変化量に重み付け係数を設定して、Ni×Ni次の対角マトリックスである重みマトリックスΠが得られる。その対角要素は、各重み付け係数からなる。ＰＭＵ有りの測定の変化量の測定品質が高いほど、また、その属するノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きい。再帰最小二乗法の計算プロセスにおける中間プロセスマトリックスの変数Pの初期値P-1はΠとし、その下付き「-1」は初期値を表すものである。ＰＭＵ有りの測定の変化量の測定品質が判明されていない場合、各重み付け係数をいずれも同じ値の０．１とする。この場合、
Figure 0005486686
（３．ｃ） 初期化忘却因子λは、０．９９とすることができる。
（４）各推定しようとする測定量Xi毎に、以下の繰り返し処理を実行する。各最新時刻（第k時間ステップ）の準定常状態断面において、対応のＳＣＡＤＡ測定値または状態推定値を利用して、以下のステップに従って係数bijを計算し、さらに、ＰＭＵ有りのノードの動態測定データを利用して、次回の準定常状態時間断面までの期間中の推定される量のダイナミクス（グラフ）を推定する。
（４．ａ）ノードiの推定しようとする測定の現在の第k時間ステップ（即ち、Tk時刻）におけるＳＣＡＤＡ計測値または状態推定値の、前回のＳＣＡＤＡまたは状態推定断面である第k-1時間ステップにおける対応値に対する変化量を算出することにより、第k時間ステップにおける出力ベクトルdkを獲得する。即ち、
Figure 0005486686
（４．ｂ）変化量ベクトルukを算出する。但し、変化量ベクトルukは、集合SPMU-i中の各ノードの、推定しようとする測定と対応する各測定の現在の第k時間ステップにおけるＳＣＡＤＡ計測値または状態推定値の、第k-1時間ステップにおける対応値に対する変化量のベクトルukである。即ち、
Figure 0005486686
（４．ｃ）次の式に従って、マトリックスgkを求める。
Figure 0005486686
（４．ｄ）次の式に従って、ベクトルwk=[ bi1(Tk) bi1(Tk) … biNi(Tk) ]Tを求める。即ち、各線形結合係数bij(Tk)を求める。
Figure 0005486686
（４．ｅ）次の式に従って、マトリックスPkを求め、次のＳＣＡＤＡまたは状態推定時刻Tk+1における線形結合係数計算のための準備を行う。
Figure 0005486686
i） 各ＰＭＵ測定時刻tにおいて、推定しようとする量XiのSPMU-i集合中の各ノードのＰＭＵ測定量の、ＳＣＡＤＡまたは状態推定の第k時間ステップの値に対する変化量であるΔXj(t)=Xj(t)−Xj(Tk) を算出し、但し、Tk<t<Tk+1。
ii）ステップ（４．ｄ）で算出したＳＣＡＤＡまたは状態推定の第k時間ステップの線形結合係数bij(Tk)を用いて、次の式
Figure 0005486686
に基づいて、ＳＣＡＤＡまたは状態推定の第k+1時間ステップになるまでの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵ無しのノードの推定しようとする量Xi(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻Tkおける該当量Xi (Tk)に対する変化量の推定値ΔXi(t)を算出する。但し、Tk<t<Tk+1 である。
iii）方程式Xi(t)=ΔXi(t)+Xi(Tk) に基づいて、当該ＰＭＵ無しのノードの推定しようとする量の、各ＰＭＵ測定時刻における推定値Xi(t)を算出し、但し、Tk<t<Tk+1。推定された相対位相角値を-180°〜+180°の範囲内に換算しなければならない。必要に応じて、推定された相対位相角値をさらに絶対位相角値に変換することができる。同様に、絶対位相角値も-180°〜+180°の範囲内に換算しなければならない。
ＳＣＡＤＡ測定または状態推定値が存在する時刻において、ＰＭＵ無しのノードの推定しようとする量は、ＳＣＡＤＡ測定値または状態推定値とする。
iv）現在のＰＭＵ測定時刻における推定しようとする量の推定値Xi(t)を記憶し、動態グラフを更新する。
v）現在のＰＭＵ測定時刻において、その他の推定しようとする量があるか否かを判断し、あると判断した場合は、ステップi）に戻り、ないと判断した場合、引き続きステップvi）を実行する。
vi）次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間、即Tk+1時刻になるまでの期間中、ＳＣＡＤＡのスイッチング量が、電力ネットワークの短路または断線故障のため、即時に更新されていることがあったか否かを検出する（通常、ＳＣＡＤＡのアナログ量は、短路または断線事故の場合に即時に更新されることがない）。上記更新があった場合、トポロジー解析を行い、ＰＭＵ無しのノードiまたはＰＭＵ有りのノードj（j∈SPMU-i）と直接接続されている線路に短路または開路が発生した場合、以下の処理を行う。
（ａ）ＰＭＵ有りのノードjが系統から孤立された場合、それと対応した線形結合係数bijをその他のＰＭＵ有りのノードに均等配分し、孤立されたＰＭＵノードjの線形結合係数を０にする。
（ｂ）ＰＭＵ無しのノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において、当該ノードのダイナミクスを推定することができないことを示す。孤立されたＰＭＵ無しのノードが再び系統に接続されたら、孤立される前の線形結合係数を初期係数bijとして、該当の測定のダイナミクスを推定する。
（ｃ）その他の状況において、即ち、ＰＭＵ無しのノードiまたはＰＭＵ有りのノードj（j∈SPMU-i）が、短路または断線のため、系統から孤立されたという事態が発生しなかった場合、特殊な処理を行わない。次回、即ち第k+1時間ステップのＳＣＡＤＡアナログ量測定または状態推定データの更新のみによって線形結合係数bijを修正して、ネットワーク構成の変化を反映する。
vii）新しい断面、即ち第k+1時間ステップのＳＣＡＤＡまたは状態推定断面が現れたか否かを判断し、現われていない場合、ステップi）に戻って、引き続き、各推定しようとする量Xiの既存の線形結合係数bij(Tk)を用いて、次のＰＭＵ測定時刻tにおける各推定しようとする量の数値を算出し、一方、現われた場合では、ステップ（４．ｆ）の繰り返しを中断して、ステップ（４．ａ）に戻って、各推定しようとする量Xiの新しい断面、即ち第k+1 時間ステップのＳＣＡＤＡまたは状態推定断面の線形結合係数bij(Tk+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を選択し、当該ノードの電圧振幅値、電圧位相角、電流振幅値、電流位相角、有効電力、無効電力、周波数および発電機回転速度を含む、当該ノードにおけるダイナミクスを推定する必要がある物理量Xiを選択する、第１のステップと、
    各ＰＭＵを備えていないノードi毎に、そのダイナミクスを推定するために使われる、ＰＭＵを備えたノードの集合SPMU-iを選択し、当該ノードの集合が、ＰＭＵを備えていないノードとの電気的距離が比較的近い複数の、ＰＭＵを備えたノードからなり、ＰＭＵを備えたノードと推定しようとするノードとが同じ電気島に属していることが保証される前提で、前記集合SPMU-i中のノードを、オフラインで形成された典型的な運転方式におけるノードインピーダンスマトリックスに基づいて選択し、または手動で指定することが可能であり、ＰＭＵを備えたノードとＰＭＵを備えていないノードとの間の最小線路段数に基づいて、オンライン自動探索方式で決定することも可能であり、集合SPMU-i中のノードの、推定しようとする測定と同じ種類の測定を、ダイナミクス推定用に選択することが可能である、第２のステップと、
    各ＰＭＵを備えていないノードiにおける推定しようとする測定量Xiの変化量ΔXiと該当のＰＭＵを備えたノードにおける対応の測定量Xjの変化量ΔXjとの間の線形結合係数bij、および再帰最小二乗法パラメータに対して、初期化を行い、但し、j∈SPMU-iとする、第３のステップと、
    各推定しようとする測定Xi毎に、現在最新のＳＣＡＤＡ測定時刻Tkにおける計測値または最新の状態推定時刻Tkにおける計算値に基づき、さらに、ヒストリー時刻におけるＳＣＡＤＡ計測値または状態推定値を参照して、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔXi(Tk)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔXj(Tk)との間の線形結合係数bij(Tk)を算出し、但し、j∈SPMU-iとし、Tkの下付き「k」が第k回のＳＣＡＤＡ測定または第k回の状態推定を示し、循環変数kが0から始まる整数を示す、第４のステップと、
    第４のステップにおいて、Tk時刻におけるＳＣＡＤＡ計測値または状態推定値に基づいて算出された、現在推定しようとする量の変化量ΔXiとＰＭＵを備えたノードにおける対応の測定量の変化量ΔXjとの間の線形結合係数bij(Tk)に基づいて、以下の式
    Figure 0005486686
     を用いて、次回のＳＣＡＤＡ測定時刻Tk+1または状態推定時刻Tk+1までの期間中の、各ＰＭＵ測定時刻tにおけるＰＭＵを備えていないノードの推定しようとする量Xi(t)の、最も近いＳＣＡＤＡ計測時刻または状態推定時刻Tkにおける該当の量Xi(Tk)に対する変化量の推定値ΔXi(t)を算出し、但し、Tk<t<Tk+1とし、さらに、方程式Xi(t)=ΔXi(t)+Xi(Tk) に基づいて、当該ＰＭＵを備えていないノードの推定しようとする量の各ＰＭＵ測定時刻における推定値Xi(t)を算出し、但し、Tk<t<Tk+1とする、第５のステップと、
    新しいＳＣＡＤＡ測定データまたは状態推定結果が得られると、第４のステップに戻り、ＰＭＵを備えていないノードの推定しようとする量の変化量とＰＭＵを備えたノードに
    おける対応の測定量の変化量との間の線形結合係数を更新する、第６のステップとを含む、ことを特徴とするダイナミクスのリアルタイム推定方法。
  2. 第２のステップにおいて、前記ＰＭＵを備えていないノードiのダイナミクスを推定するために使われる、ＰＭＵを備えたノードの集合SPMU-iは、選択された４〜６つのノードから構成されている、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  3. 第２のステップにおいて、前記オンライン自動探索方法は、トポロジー情報に基づいて、各ＰＭＵノードと推定しようとするノードiとの間の最小線路段数を自動的に計算し、ノードiと１段の線路を介して接続、即ち直接接続しているＰＭＵノードを集合SPMU-iの中に優先に入れてから、前記集合中のノード数が予め決定された数になるまで、ノードiとそれぞれ２段の線路、３段の線路、４段の線路を介して接続しているＰＭＵを備えたノードを、介した線路の段数が小さい方から順に集合SPMU-iに加える方法である、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  4. 第２のステップにおいて、推定しようとする測定量が電圧振幅値または電流振幅値である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、両電圧レベルの定格電圧比を用いて、ＰＭＵを備えたノードの電圧振幅値または電流振幅値を、ＰＭＵを備えていないノードの電圧レベルにおける電圧振幅値または電流振幅値に変換する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  5. 第２のステップにおいて、推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵによって直接測定した絶対位相角を、当該電気島に存在するノードの電圧フェーザを基準とした相対位相角に変換し、即ち、相対位相角に基づいてダイナミクスの推定を行い、かつ、相対位相角を-180°〜+180°の範囲内に換算する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  6. 第２のステップにおいて、推定しようとする測定量が電圧位相角または電流位相角である場合、ＰＭＵを備えたノードとＰＭＵを備えていないノードとは電圧レベルが異なっているならば、さらに、変圧器の結線方式によって生じた位相シフトに基づいて、ＰＭＵを備えたノードの電圧位相角または電流位相角を、ＰＭＵを備えていないノードの電圧レベルのものに変換し、かつ、-180°〜+180°の範囲内に換算する必要がある、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  7. 第３のステップにおける初期化の内容は、下記の処理（ａ）、処理（ｂ）、処理（ｃ）、処理（ｄ）をさらに含み、
    処理（ａ）：まず、推定しようとする測定量の変化量ΔXiの線形結合係数bijを求めるという問題の再帰最小二乗法の数学モデルを作成し、再帰最小二乗法を用いて問題を解く場合の一般的な形式は、
    Figure 0005486686
     であり、dが観測された出力ベクトル、wが求めようとするパラメータベクトル、Aが複数のベクトルujからなる状態マトリックスであり、また、ujがパラメータベクトルwと観測量d( j )との関係を表し、vがノイズベクトルであり、推定しようとする測定量の変化量ΔXiの線形結合係数bijを求めるという問題について、前記式において、計算中に必要な変数の定義は、
    Figure 0005486686
     となっており、ΔXiは、推定しようとする測定Xiの現在のＳＣＡＤＡ測定値と前回のＳＣＡＤＡ測定値または状態推定値との差値、または推定しようとする測定Xiの現在の状態推
    定値と前回の状態推定値との差値であり、各推定しようとする測定をそれぞれ推定するため、実際には、出力ベクトルdは１つの要素しか存在せず、ΔX1 ΔX2 … ΔXj … ΔXNiは、それぞれ、SPMU-i中の各ノードの現在のＳＣＡＤＡ測定値または状態推定値と前回のＳＣＡＤＡ測定値または状態推定値との差値であり、bi1 bi2 … bij … biNi は、それぞれ、SPMU-i中の各ノードの該当の測定量の変化量とノードiにおける推定しようとする測定量の変化量ΔXiとの間の線形結合係数であり、
    処理（ｂ）：前記推定しようとする測定量の変化量ΔXiの線形結合係数 bijの初期値を1/Niとすると、再帰最小二乗法で求めようとするパラメータベクトルwの初期値はw-1=[1/Ni 1/Ni … 1/Ni ]Tとなり、但し、NiはＰＭＵを備えていないノードiのSPMU-i集合中の要素の数であり、下付き「-1」は初期値を表すものであり、
    処理（ｃ）：前記各線形結合係数と対応する、ＰＭＵを備えた測定の変化量に対して、重み付け係数を設定することにより、Ni×Ni次の対角マトリックスとなる重みマトリックスΠが得られ、その対角要素は各重み付け係数からなり、ＰＭＵを備えた測定の変化量の測定品質が高いほど、また、そのノードと推定されるノードとの電気的距離が近いほど、対応の重み付け係数が大きく、再帰最小二乗法の計算プロセスにおける中間プロセスマトリックス変数Pの初期値P-1をΠとし、ＰＭＵを備えた測定の変化量の測定品質が不明な場合、各重み付け係数を同じ値の０．１とし、この場合、
    Figure 0005486686
     となり、但し、下付「-1」は初期値を表すものであり、
    処理（ｄ）： 望ましくは、初期化忘却因子λを０．９９とする、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  8. 第４のステップにおいて、再帰最小二乗法を用いて、ＰＭＵを備えていないノードの推定しようとする量の変化量ΔXi(Tk)とＰＭＵを備えたノードにおける対応の測定量の変化量ΔXj(Tk)との間の線形結合係数bij(Tk)を求める計算には、下記のステップ（ａ）、ステップ（ｂ）、ステップ（ｃ）を含み、
    ステップ（ａ）：ノードiの推定しようとする測定の現在の第k時間ステップ、即ちTk時刻におけるＳＣＡＤＡ測定値または状態推定値の、前回のサンプリング時刻におけるＳＣＡＤＡ測定値または状態推定値である第k-1時間ステップにおける対応値に対する変化量を算出することにより、第k時間ステップの出力ベクトルdkが得られ、即ち、
    Figure 0005486686
     ステップ（ｂ）：次の式
    Figure 0005486686
     により、前記集合S PMU-i 中の各ノードの、推定しようとする測定と対応する各測定の現在の第k時間ステップにおけるＳＣＡＤＡ測定値または状態推定値の、第k-1時間ステップにおける対応値に対する変化量のベクトルu k を算出し、
    ステップ（ｃ）：次の３つの回帰公式
    Figure 0005486686
     に基づいて、第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時間ステップの各線形組み合わせ係数bij(Tk)が得られる、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  9. 第５のステップにおいて、次回のＳＣＡＤＡ測定時刻または状態推定時刻を待っている間に、即時に更新されたＳＣＡＤＡスイッチング量が、電力ネットワークに短絡事故または断線事故が発生したことを示した場合、さらに、ＰＭＵを備えていないノードiまたはＰＭＵを備えたノードj（j∈SPMU-i）に直接接続されている線路に、短路または断線が発生した場合、下記の処理（ａ）、処理（ｂ）および処理（ｃ）を実施し、
    処理（ａ）：ＰＭＵを備えたノードjが系統から孤立された場合、それと対応した線形結合係数bijをその他のＰＭＵを備えたノードに平均配分し、孤立されたＰＭＵを備えたノードjの線形結合係数bijを0とし、
    処理（ｂ）：ＰＭＵを備えていないノードiが系統から孤立された場合、それと対応した測定値を０にすることにより、孤立期間において当該ノードのダイナミクスを推定することが不可能であることを示し、孤立されたＰＭＵを備えていないノードiが再び系統に接続されたら、孤立される前の線形結合係数を初期係数bijとして、該当の測定のダイナミクスを推定し、
    処理（ｃ）：その他の状況、即ち、短絡または断線によってＰＭＵを備えていないノードiまたはＰＭＵを備えたノードjが系統から孤立されたという状況にならなかった場合、特殊な処理を行わず、次回の、即ち第k+1時間ステップのＳＣＡＤＡアナログ量測定データまたは状態推定データの更新によって、線形結合係数bijを補正して、ネットワーク構成の変化を表す、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
  10. 第５のステップにおいて、推定された相対位相角値を-180°〜+180°の範囲内に換算し、または、推定された相対位相角値を絶対位相角値に変換し、さらに、-180°〜+180°の範囲内に換算する、ことを特徴とする請求項１に記載のダイナミクスのリアルタイム推定方法。
JP2012530072A 2009-09-25 2009-12-09 状態マトリックスに依存せずにｐｍｕ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法 Active JP5486686B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
CN2009100930067A CN101661069B (zh) 2009-09-25 2009-09-25 不依赖状态矩阵的弱可观非pmu测点动态过程实时估计方法
CN200910093006.7 2009-09-25
PCT/CN2009/001407 WO2011035461A1 (zh) 2009-09-25 2009-12-09 不依赖状态矩阵的弱可观非pmu测点动态过程实时估计方法

Publications (2)

Publication Number Publication Date
JP2013506389A JP2013506389A (ja) 2013-02-21
JP5486686B2 true JP5486686B2 (ja) 2014-05-07

Family

ID=41789207

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012530072A Active JP5486686B2 (ja) 2009-09-25 2009-12-09 状態マトリックスに依存せずにｐｍｕ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法

Country Status (5)

Country Link
US (1) US9178386B2 (ja)
JP (1) JP5486686B2 (ja)
CN (1) CN101661069B (ja)
AU (1) AU2009353270B2 (ja)
WO (1) WO2011035461A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US12100958B2 (en) 2019-08-29 2024-09-24 Hitachi Energy Ltd Cloud-based hybrid state estimation

Families Citing this family (58)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2580696A2 (en) * 2010-06-11 2013-04-17 ABB Research Ltd. Detecting state estimation network model data errors
CN102361353B (zh) * 2011-10-26 2013-09-25 东北电网有限公司 基于双时间尺度延时估计的断面生数据对齐方法
US9627886B2 (en) * 2012-03-27 2017-04-18 Mitsubishi Electric Research Laboratoriies, Inc. State estimation for power system using hybrid measurements
US9337662B2 (en) 2012-06-15 2016-05-10 Abb Technology Ag Parallel computation of dynamic state estimation for power system
US9989580B2 (en) * 2012-08-14 2018-06-05 Abb Research Ltd. Fault analysis in electric networks having a plurality of multi-phase buses
CN102832617B (zh) * 2012-09-03 2014-08-13 武汉大学 一种基于精确化模式判别的大电网暂态稳定分析方法
CN103076556B (zh) * 2013-01-21 2013-10-16 中国民航大学 航电组件功能维护测试点选取方法
CN103116097B (zh) * 2013-01-25 2015-06-24 中国电力科学研究院 基于多断面混合量测信息的设备参数在线辨识方法
US9488682B2 (en) * 2013-02-20 2016-11-08 Mitsubishi Electric Research Laboratories, Inc. Method for determining power consumption of loads in ungrounded power distribution systems
US9876352B2 (en) * 2013-02-26 2018-01-23 Washington State University Voltage stability monitoring in power systems
CN103198437B (zh) * 2013-04-18 2016-05-11 国家电网公司 一种电网量测数据和电网模型纠错方法及装置
CN103326358B (zh) * 2013-06-17 2015-04-29 西南交通大学 基于同步相角测量装置的电力系统动态状态估计方法
CN103390251B (zh) * 2013-07-24 2016-06-01 大连海事大学 一种电力系统状态估计的量测量权重设置方法
CN103413044B (zh) * 2013-08-13 2016-03-30 东北电网有限公司 一种基于变电站量测信息的电力系统局部拓扑估计方法
WO2015052747A1 (ja) * 2013-10-07 2015-04-16 株式会社日立製作所 電力系統状態推定装置およびその電力系統状態推定方法
US20150168465A1 (en) * 2013-12-13 2015-06-18 General Electric Company Method and apparatus for electric power system distribution state estimations
CN103838962B (zh) * 2014-02-18 2017-02-22 河海大学 一种计及pmu量测的分步线性状态估计方法
CN103902814B (zh) * 2014-03-10 2017-06-16 中国南方电网有限责任公司 基于动态分区的电力系统运行状态检测方法
US20150346703A1 (en) * 2014-05-27 2015-12-03 Infineon Technologies Ag State observers
EP3032270A1 (en) * 2014-12-12 2016-06-15 Airbus Defence and Space, S.A. Method and system for performing electrical tests to complex devices
CN105046038A (zh) * 2015-04-03 2015-11-11 国家电网公司 一种线路电纳参数的估计方法
CN105139286A (zh) * 2015-08-31 2015-12-09 国网山东省电力公司菏泽供电公司 基于无功不平衡量的scada数据有效性评估方法
CN105759169B (zh) * 2016-03-07 2018-07-06 华中科技大学 一种基于wams时间断面信息和拓扑信息的故障诊断方法
CN106099906B (zh) * 2016-05-13 2018-09-14 国电南瑞科技股份有限公司 一种电力系统安全稳定在线综合防御方法
CN106053930B (zh) * 2016-07-06 2018-11-23 西北核技术研究所 一种抗随机噪声的无开关弱信号同步积分测量装置及测量方法
CN106208041B (zh) * 2016-07-18 2019-03-01 西南交通大学 一种基于分段常数独立分量分析的配电网谐波电流估计方法
CN106405267B (zh) * 2016-07-27 2018-10-19 河北省电力建设调整试验所 一种广域后备保护系统中的故障检测方法
CN106296461B (zh) * 2016-08-17 2020-05-12 华北电力大学 基于粒子群优化算法及局部状态估计的电网参数辨识方法
CN106355029B (zh) * 2016-09-14 2019-09-27 国网天津市电力公司 一种配电系统关键量测解耦快速辨识的方法
CN106646106B (zh) * 2016-10-11 2019-02-22 河海大学 基于变点探测技术的电网故障检测方法
CN106570345A (zh) * 2016-11-15 2017-04-19 中国电力科学研究院 一种基于图论理论的参数辨识方法及装置
CN107046285A (zh) * 2017-04-12 2017-08-15 国家电网公司 一种基于混合量测的配电网状态评估方法
CN107271768B (zh) * 2017-05-26 2019-06-21 东南大学 一种最小二乘拟合动态频率测量方法
CN109217282B (zh) * 2017-06-29 2024-07-02 中国电力工程顾问集团华东电力设计院有限公司 一种确定主变环流运行参数的系统及方法
CN107425523B (zh) * 2017-08-09 2020-05-22 国网辽宁省电力有限公司 一种复杂电力系统运行误差识别与自动校正方法及系统
JP6858271B2 (ja) * 2017-10-31 2021-04-14 株式会社日立製作所 電力系統の安定性解析装置、安定化装置および方法
CN109031000B (zh) * 2018-08-03 2019-08-16 贵州电网有限责任公司电网规划研究中心 一种基于非故障扰动就地测量电网短路容量的方法及系统
CN109472295B (zh) * 2018-10-17 2021-11-12 贵州电网有限责任公司 基于Pearson相关系数的WAMS/SCADA系统对时排序方法
CN109921415B (zh) * 2019-03-13 2022-08-02 东北大学 一种面向混合量测的电网恶性数据注入攻击在线防御方法
CN109993361B (zh) * 2019-03-29 2023-06-13 国网上海市电力公司 一种基于pmu的配电网运行趋势预测方法
CN110210690B (zh) * 2019-06-21 2022-12-06 国网天津市电力公司 一种配电系统微型同步相量测量单元优化配置方法
CN110378041B (zh) * 2019-07-23 2022-05-27 广东电网有限责任公司 一种电网状态估计方法、系统、设备及计算机存储介质
US11119462B2 (en) * 2019-10-08 2021-09-14 Global Energy Interconnection Research Institute (Gelrina) Systems and methods for hybrid dynamic state estimation
CN110867854B (zh) * 2019-11-26 2021-10-08 广东电网有限责任公司广州供电局 配网可观测性快速检测方法
CN111551820B (zh) * 2020-04-21 2022-03-08 天津大学 基于同步相量量测的配电网快速频率调节的系统辨识方法
CN112003317B (zh) * 2020-06-05 2022-01-25 国电南瑞科技股份有限公司 一种水电机组功率调节输入输出优化方法及装置
CN111600393B (zh) * 2020-06-18 2020-12-25 国网四川省电力公司电力科学研究院 变电站不同电压等级的电压量测数据归算方法及识别方法
CN111884346B (zh) * 2020-08-17 2021-10-22 国网四川省电力公司电力科学研究院 基于广域量测信息的低压配电网保护方法和介质
CN112115568B (zh) * 2020-09-01 2023-12-01 贵州电网有限责任公司 一种基于pmu量测的改进耦合最小二乘输电线路参数辨识方法
CN112115593B (zh) * 2020-09-03 2023-03-10 西北工业大学 一种基于中心对称多面体的分布式集员估计方法
CN112886579B (zh) * 2021-02-03 2022-03-25 北京交通大学 考虑pmu和scada量测时标不匹配的配电网线路参数辨识方法
CN112837180B (zh) * 2021-02-05 2022-06-03 国网湖南省电力有限公司 用于电力系统的scada-pmu数据融合方法
CN113095564A (zh) * 2021-04-08 2021-07-09 国网河北省电力有限公司经济技术研究院 一种配电网预测辅助装置及其预测状态估计方法
CN114002512A (zh) * 2021-10-29 2022-02-01 华中科技大学 一种机电瞬时频率动态的测量计算方法及装置
CN115051382B (zh) * 2022-06-22 2024-03-29 西南交通大学 机车接入三相电网低频振荡模态分析导纳模型的建立方法
CN117559454B (zh) * 2023-11-10 2024-05-31 哈尔滨工业大学 一种信息-物理模型联合驱动的电网拓扑结构辨识方法
CN117526318B (zh) * 2024-01-05 2024-05-10 珠海许继电气有限公司 配电网状态评估方法、装置、电子设备、存储介质
CN118336929B (zh) * 2024-06-13 2024-08-20 广东电网有限责任公司梅州供电局 一种电网状态监控方法、装置、系统和介质

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE60108966D1 (de) * 2001-12-21 2005-03-24 Abb Research Ltd Zustandsbestimmung für ein Energieverteilungsnetz
DE60119555T2 (de) * 2001-12-21 2007-03-08 Abb Schweiz Ag Ermittlung der Betriebsgrenzwerte in einem Energieverteilungsnetz
ES2299678T3 (es) * 2003-06-21 2008-06-01 Abb Research Ltd. Metodo para detectar oscilaciones electromecanicas en sistema de potencia.
US8190299B2 (en) * 2006-07-19 2012-05-29 Rovnyak Steven M Integrated and optimized distributed generation and interconnect system controller
JP2008154362A (ja) * 2006-12-18 2008-07-03 Tokyo Electric Power Co Inc:The 電力系統の状態推定装置及び方法
CN100486073C (zh) * 2007-05-15 2009-05-06 浙江大学 考虑负荷电压特性的非线性动态状态估计算法
CN101291061B (zh) * 2008-05-16 2010-04-21 南京南瑞继保电气有限公司 一种电力系统动态过程状态估计方法
US7979220B2 (en) * 2008-08-21 2011-07-12 Abb Research Ltd. System and method for monitoring power damping compliance of a power generation unit
CN101383511B (zh) * 2008-10-10 2010-08-04 清华大学 基于数据采集系统量测数据的电力系统状态估计方法
CN101505059B (zh) * 2009-03-06 2011-01-26 南京南瑞继保电气有限公司 调度能量管理系统中输电网电抗参数估计方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US12100958B2 (en) 2019-08-29 2024-09-24 Hitachi Energy Ltd Cloud-based hybrid state estimation

Also Published As

Publication number Publication date
CN101661069A (zh) 2010-03-03
AU2009353270B2 (en) 2015-04-30
CN101661069B (zh) 2011-07-20
WO2011035461A1 (zh) 2011-03-31
JP2013506389A (ja) 2013-02-21
AU2009353270A1 (en) 2012-03-08
US20120283967A1 (en) 2012-11-08
US9178386B2 (en) 2015-11-03

Similar Documents

Publication Publication Date Title
JP5486686B2 (ja) 状態マトリックスに依存せずにｐｍｕ無しの低可観測性測定点のダイナミクスをリアルタイムで推定する方法
CN107843810B (zh) 一种基于状态估计的有源配电网故障区段在线定位方法
Schenato et al. Bayesian linear state estimation using smart meters and pmus measurements in distribution grids
Della Giustina et al. Electrical distribution system state estimation: measurement issues and challenges
EP1946125B1 (en) Method and apparatus for verifying the accuracy of state estimation calculations
CN110299762B (zh) 基于pmu准实时数据的主动配电网抗差估计方法
Zhao et al. On-line PMU-based transmission line parameter identification
CN103020726B (zh) 面向全pmu量测的抗差状态估计方法
US11016134B2 (en) Determination of power transmission line parameters using asynchronous measurements
Atanackovic et al. Deployment of real-time state estimator and load flow in BC Hydro DMS-challenges and opportunities
CN103983847B (zh) 一种同步相量测量中基于rls的自适应频率跟踪测量方法
CN109754013A (zh) 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法
CN102522742B (zh) 基于单点量测信息的外网戴维南等值参数的估计方法
Watitwa et al. A review on active distribution system state estimation
CN114065118B (zh) 一种基于指数函数的电力系统鲁棒状态估计方法
CN109494733B (zh) 一种电力负荷模型的辨识参数优化方法及系统
CN109858061B (zh) 面向电压功率灵敏度估计的配电网络等值与化简方法
CN108896852B (zh) 一种公共接入点短路容量在线测量方法及系统
Angioni et al. Impact of heterogeneous measurements in the state estimation of unbalanced distribution networks
Goklani et al. Quantification of minimum unbalance required for accurate estimation of sequence parameters of transmission line using PMU data
Baran et al. State estimation for real time monitoring of distribution feeders
Mingorança et al. Critical buses identification for voltage stability assessment considering the application of modal analysis and a robust state estimation with bad data suppression
CN108736475B (zh) 基于pmu监测的互联电网子系统运行可靠性评估方法
Lu et al. PMU assisted power system parameter calibration at Jiangsu electric power company
Liu et al. A New Method of Data Fusion with PMU and SCADA Based on Branch Decomposition

Legal Events

Date Code Title Description
A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20131119

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20140107

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20140218

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20140221

R150 Certificate of patent or registration of utility model

Ref document number: 5486686

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250