JPWO2015136666A1 - 個別電気機器稼働状態推定装置、およびその方法 - Google Patents

Abstract

個別電気機器稼働状態推定装置は、総使用電力測定手段（5）と、総使用電力の時間的変動（ジャンプ電力）を算出する時間変動算出手段（6a）と、ジャンプ電力から新規電気機器として既存の機種リストに加えるか、また余分な機種を削除するかを行う電気機器機種特定手段（6b）と、アップおよびダウのジャンプがどの電気機器のオン/オフの事象の生起に対応しているかの尤度を推定する対応尤度推定手段（6c）と、尤度が入力されたときの個別電気機器の稼働状態の変化を推定して現在の稼働状態の推定を更新し、これと規格電力の推定結果を用いて総使用電力の予測値と総使用電力の測定データから個別電気機器の機種、規格電力、および稼働状態の各推定を最適化し、個別電気機器の稼働確率を動的に推定する個別電気機器稼働状態推定手段（6d）と、を有する。

Description

本発明は、複数の電気機器のうち、稼働している電気機器の種類と、その稼働状態とを同時にリアルタイムでモニタリング可能な個別電気機器稼働状態推定装置、およびその方法に関する。

電力の効率的な利用を図るため、電気機器について個別にそれらの使用状況を把握できるようにすることが望ましい。この場合、各電気機器にセンサを取付ければ、個々の電気機器の使用状況（電気機器ごとの稼働・非稼働の状態やその電力使用量）が個別に測定でき、この測定結果に基づいて電気機器の使用状況を把握することで、省エネルギ化や電力供給制御等に役立てることが可能となる。
ところが、このように各電気機器にセンサをそれぞれ設けると、センサ数が増えるのみならず、それらの取付けコストや、センサからの測定信号の収集システム（無線等）構築のためのコストが必要となり、高コスト化をまぬがれない。
そこで、安価な装置として一般家庭等などでも利用可能とするため、各電気機器にそれぞれセンサを設けなくとも各電気機器の使用状況をモニタリング可能とした技術が開発されてきている。

このような従来の個別電気機器稼働状態推定装置およびその方法は、特許文献１に記載の技術が知られている。
この特許文献１に記載の個別電気機器稼働状態推定装置では、使用している電気機器の総使用電力の電気特性に基づきパルス型波形のパターン認識を利用して、電気機器の電力使用や操作特性を個別に測定する非侵入型電力使用モニタリング（NIALM: Non-intrusive Appliance Load Monitoring）を行うようにしている。

すなわち、より具体的には、その従来の個別電気機器稼働状態推定装置は、設備の回路での全負荷の電気特性（コンダクタンスやサセプタンスからなるアドミッタンス）を測定し、それらに比例したアナログ信号を生成し、設備外の位置で上記回路に接続されたセンサ手段（電圧センサ、電流センサ）と、測定した電圧や電流からコンダクタンスやサセプタンスを演算し、上記アナログ信号をディジタル信号に変換するＡＤ変換器と、ディジタル信号に応じて複数のユニットの各々による全負荷をそれらの部分に分ける信号プロセッサとを備えている。なお、この処理結果を利用して、個別に電気機器の稼働状況をディスプレイに表示可能としている。

米国特許第４，８５８，１４１号公報

しかしながら、上記従来の個別電気機器稼働状態推定装置には、以下に説明するような問題がある。
上記従来技術にあっては、その適用にあたって、どの家庭にどのような規格電力を持ったどのような電気機器が何台保有されているかといった情報を事前に収集したり、また廃棄される電気機器や新たに追加される電気機器がある場合の情報をその都度更新したりする必要がある。このため、大きな手間やコストがかかる上、一般家庭で素人がそれらの情報を入力したり更新したりすることは期待できず、現実には適用が非常に難しいといった問題がある。

本発明は、上記問題に着目してなされたもので、その目的とするところは、個々の電気機器にセンサを設けることなく、また事前にどのような種類の電気機器が何台あるかといった情報がなくても、電気機器の種類の特定とそれらの稼働状態とをリアルタイムでモニタリングすることができるようにした個別電気機器稼働状態推定装置、およびその方法を提供することにある。

この目的を達成するため、本発明による個別電気機器稼働状態推定装置およびその方法では、統計的逆推定の方法によって、総使用電力からリアルタイムに個別電気機器の機種とその稼働状態とを推定する動的アルゴリズムを構築する。このため、総使用電力の時間的変動から、発生頻度の高い変動(ジャンプ電力)をクラスタリングして個別電気機器（機種）ごとの規格電力を推定し、機種のクラスタと計測されたジャンプとからリアルタイムで電気機器の特定およびその稼働状況を推定するとともに、必要に応じてクラスタの変更を行うアルゴリズムを実行するようにしている。

すなわち、本発明の個別電気機器稼働状態推定装置は、
複数の電気機器の総使用電力を測定する総使用電力測定手段と、
総使用電力測定手段で測定した総使用電力に基づき、この総使用電力の時間的変動であるジャンプ電力を算出する時間変動算出手段と、
時間変動算出手段で算出したジャンプ電力に基づいて新たな電気機器として既存の機種リストに加えるか否かを判定し、新たな電気機器として既存の機種リストに加える場合には個別電気機器別の規格電力の再計算を行い、必要に応じて既存の機種リストから余分な機種を削除し、新たな電気機器を加えない場合には既存の機種リストを維持するクラスタリングを行う電気機器機種特定手段と、
時間変動算出手段で算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定する対応尤度推定手段と、
対応尤度推定手段で得られた、現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの個別電気機器の稼働状態の変化を推定して、現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と総使用電力の測定データから個別電気機器の機種の推定、規格電力の推定、および稼働状態の推定を最適化し、測定された総使用電力の時間変動の時間間隔で刻々変化する個別電気機器の稼働確率を動的に推定する個別電気機器稼働状態推定手段と、
を備えたことを特徴とする。

一方、本発明の個別電気機器稼働状態推定方法は、
複数の電気機器の総使用電力を測定するステップと、
測定した総使用電力に基づき、この総使用電力の時間的変動であるジャンプ電力を算出するステップと、
算出したジャンプ電力に基づいて新たな電気機器として既存の機種リストに加えるか否かを判定し、新たな電気機器として既存の機種リストに加える場合には個別電気機器別の規格電力の再計算を行い、必要に応じて既存の機種リストから余分な機種を削除し、新たな電気機器を加えない場合には既存の機種リストを維持するクラスタリングを行うステップと、
算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定するステップと、
現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの個別電気機器の稼働状態の変化を推定して、現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と総使用電力の測定データから個別電気機器の機種の推定、規格電力の推定、および稼働状態の推定を最適化し、測定された総使用電力の時間変動の時間間隔で刻々変化する個別電気機器の稼働確率を動的に推定するステップと、
を有することを特徴とする。

本発明の個別電気機器稼働状態推定装置およびその方法にあっては、個々の電気機器にセンサを設けることなく、また事前にどのような種類の電気機器が何台あるかといった情報がなくても、モニタリングの対象となる個々の電気機器の使用状況をリアルタイムでモニタすることができる。

複数の電気機器を備えた家屋内に、本発明の実施例１の個別電気機器稼働状態推定装置を設けた状態を示す図である。 実施例１の個別電気機器稼働状態推定装置の構成を示す機能ブロック図である。 実施例１の個別電気機器稼働状態推定装置で実行する個別電気機器稼働状態推定の流れの考え方を表したフローチャートを示す図である。 実施例１の個別電気機器稼働状態推定装置で実行するk-ミーンズアルゴリズムを説明する（k＝2の場合）ための模式図である。 実施例１の個別電気機器稼働状態推定装置で用いるマルコフスイッチングモデルでジャンプが正の場合の状態変化をまとめた表を示す図である。 実施例１の個別電気機器稼働状態推定装置で実行するマルコフスイッチングモデルでジャンプが負の場合の状態変化をまとめた表を示す図である。 実施例１の個別電気機器稼働状態推定装置で実行する個別電気機器稼働状態推定の流れを数式とともに表したフローチャートを示す図である。 実施例１の個別電気機器稼働状態推定装置の効果を確認する実験での、実際の家庭環境データの一部としての天気のデータを示す図である。 上記実験での、総使用電力の１分間のジャンプ量の規模別度数分布を示す図である。 上記実験での、総使用電力規模別の総電力使用量に占めるシェアを示す図である。 上記実験での、８日間実環境データのリアルタイム処理アルゴリズムの結果（使用開始から11520分まで）を示す図である。 図１１の結果のうち使用開始から150分までの部分を拡大した図である。 図１１の結果のうち150分から350分までの部分を拡大した図である。 個別電気機器を手動でオンとオフした疑似家庭環境での４時間の１分間隔のデータを用いて、本発明を適用し、個別電気機器の種類および規格電力を推定し、総使用電力と個別電気機器の稼働状態を推定した結果を示す図である。 個別電気機器の種類および規格電力を推定し、個別電気機器の種類の数が変わったときの状態変数の引き渡しを考慮した上で、総使用電力と個別電気機器の稼働状態を推定した結果を示す図である。

まず、本発明の個別電気機器稼働状態推定装置およびその方法は、下記の点で、上記特許文献１に記載の技術と異なる。すなわち、これらは、
(a)総使用電力の時間変動であるジャンプ、このうち特に各電気機器のオン-オフのみに着目しているか、
(b) k-ミーンズ法（k -Means Clustering）のアルゴリズムなどの統計的アプローチをとっているか、
(c)どのような電気機器が稼働しているのかなどの事前情報をどの程度必要としているか、
(d)１分間隔の計測など、実際の環境での総使用電力の計測時系列データに対応できる動的なモデル化が確立されているか、
といった点で異なっている。
これらの性質をすべて備えたアルゴリズムとして本発明では逆推定アルゴリズムを用いるが、現時点では上記すべての性質を満たすNIALMの解法アルゴリズムは存在していない。
以下、本発明の実施の形態を、図面に示す実施例に基づき詳細に説明する。

実施例１では、本実施例の個別電気機器稼働推定装置を一般家屋(家庭)に適用した例について説明する。
図１は、上記のような逆推定アルゴリズムを用いた本実施例１の個別電気機器稼働状態推定装置が、このモニタリングの対象となる複数の電気機器を備えた家屋に設置された状態の全体構成を模式的に示す。なお、逆推定アルゴリズムとは、不完全なデータのもとで完全データを統計的に復元する方法である。
同図において、家屋１内には、複数の配線2が分電盤8から張り巡らされており、これらの配線2の適宜箇所には複数機種の電気機器3a〜3n（nは、２以上の正の整数）がそれぞれタップやコンセント等を介して接続、あるいは接続可能である。このような電気機器3a〜3nとしては、たとえば、冷蔵庫、テレビ、シーリングライト、エアコン(分電盤8を介して前者等とは独立系統で接続)などのようにオン/オフ可能な家庭用電気機器の他に、ウォッシュレット（商品名）のような不定期で負荷が発生する電気機器や、インターフォンなどのように常時オンとされている電気機器が含まれる。

家屋1内の分電盤8と屋外の供給電線4との間には、電力計5が設けられており、その家屋1で使用している総使用電力を測定する。また、分電盤8にてブレーカごとの電力を測定するものもある。これら電力計5や分電盤8は、本発明の総使用電力測定手段に相当する。以下は、個別電気機器稼働状態推定装置が電力計5で測定した総使用電力を利用する場合について説明するが、分電盤8でその分電盤ごとの総使用電力を測定するようにしてもよい。
この電力計5で計測した総使用電力に関する信号は、稼働推定部6に送られ、ここでどの電気機器が稼働しているのかといった電気機器の機種の特定、またその稼働状態について推測し、その推測結果をディスプレイ7に表示するように構成している。電力計5、稼働推定部6、およびディスプレイ7は、本実施例の個別電気機器稼働推定装置を構成する。

上記稼働推定部6は、たとえばマイクロコンピュータで構成され、図２にその構成を示す。
すなわち、図２に示すように、稼働推定部6は、時間変動算出部6aと、電気機器機種特定部6bと、対応尤度推定部6cと、個別電気機器稼働状態推定部6dと、電気機器のリスト6gと、を備えている。
以下、これらについてそれぞれ説明する。

時間変動算出部6aは、電力計5で計測した、家屋1での総使用電力に関する信号が入力されて、この信号を所定時間（本実施例では１分）間隔で現在とそのひとつ前の総使用電力の時間変動量、すなわち総使用電力の時間変動１段階差を算出し、稼働電気機器機種特定部6bと、対応尤度推定部6cと、個別電気機器稼働状態推定部6dとへ出力する。
なお、時間変動算出部6aは、本発明の時間変動算出手段に相当する。

電気機器機種特定部6bは、時間変動算出部6aで算出した時間変動1段階差（ジャンプに相当）に基づいて個別電気機器の規格電力を推定し、発生頻度の高い変動量（ジャンプ電力）をクラスタリングし、ジャンプ電力のクラスタ(機種)ごとの平均をそれぞれの個別電気器機器の規格電力として推定する。
ここで、新しい機種の判定基準では、ジャンプ量が、クラスタに対応する規格電力の分散のある一定比率を超えた場合に新しい機種として判定し、そうでない場合には既存の機種の電力使用が規格電力の通常の変動範囲内に入っていると判定する。

すなわち、電気機器機種特定部6bは、新しい電気機器があるか否かを抽出し、新しい電気機器が抽出された場合には新たな電気機器の機種として既存の機種リスト6gに加え、新しい機種が抽出されないときは既存の機種リストを維持する。
また、新しい電気機器が抽出された場合には、個別電気機器の機種別の規格電力を再計算し、必要に応じて、実質重複する機種など、余分な機種を機種リスト6gから削除（枝刈り）する。ここで、機種のクラスタ数は、可変としてある。ジャンプに応じた個別電気機器の規格電力の推定は、リアルタイムで学習していく。このようにして得た新しい機種リスト6gのデータを、対応尤度推定部6cと、個別電気機器稼働状態推定部6dと、機種リスト6gとへそれぞれ出力する。なお、電気機器機種特定部6bは、本発明の電気機器機種特定手段に相当する。

対応尤度推定部6cは、時間変動算出部6aで算出した時間変動1段階差、すなわちアップ（計測したジャンプの符号が正）とダウン（計測したジャンプの符号が負）とのジャンプが、どの電気機器のスイッチオンとスイッチオフの事象の生起に対応しているかの尤度（確率）を、電気機器機種特定部6bで得た最新の機種リストの機種ごとに推定する。この尤度は、個別電気機器稼働状態推定部6dへ出力する。
なお、対応尤度推定部6cは、本発明の対応尤度推定手段に相当する。

個別電気機器稼働状態推定部6dは、時間変動算出部6aで算出した時間変動1段階差（ジャンプ）と対応尤度推定部6cで推定した尤度とに基づいて、どの電気機器が稼働しているか、およびその稼働状況をリアルタイムで学習していく。そのため、個別電気機器稼働状態推定部6dは、マルコフスイッチングモデル6eとカルマンフィルタ6fとを備えている。
すなわち、前者を用いて現在の個別電気機器の稼働状態のもので総使用電力の変動ジャンプが得られる個別電気機器のスイッチオン／スイッチオフの事象生成に関する尤度（確率）が入力されたとき、どのような個別電気機器の稼働状態が変化するか、を推定し、後者を用いて総使用電力の時間変動の測定データと上記稼働状態の変化とから、個別電気機器の機種推定、規格電力推定および稼働状態の推定を最適化する。そして、個別電気機器の稼働状態の推定結果を、平日・土日別、時間帯別に集計していき、個別電気機器の使用パターンの学習を行う。個別電気機器の稼働状態の推定結果およびその過去のパターンから事前分布を改定していくため、この学習結果を対応尤度推定部6cやディスプレイ7へ出力する。
なお、個別電気機器稼働状態推定部6dは、本発明の個別電気機器稼働状態推定手段に相当する。

電気機器のリスト6gは、抽出された機種の異なった規格電力の値（平均値、分散）を記憶する。ここで記憶されたデータは、電気機器機種特定部6bに送られてここでの演算に用いられ、演算の結果を受けて必要に応じてデータを書き換える。

次に、上記稼働状態推定部6で実行される個別電気機器の稼働状態推定プロセスの流れについて、図３のフローチャートに基づいて説明する。
個別電気機器の稼働状態推定にあっては、上述のように、NIALMリアルタイム逆推定アルゴリズムを用いる。ここで、総使用電力の時間変動を表現している動的推定モデルの核となる部分は、一般的に、状態空間モデル(State Space Model)と呼ばれているモデルである。
この逆推定アルゴリズムの最終的なアウトプットは、測定された総使用電力の時間変動の入力ストリームに対応した時間間隔(本実施例では1分)で、刻々変化する個別電気機器のスイッチオンとスイッチオフの稼働状態を推定することである。
このように総使用電力の時間変動としてのジャンプに着目しているのは、電気機器がスイッチオンしたときに規格電力の垂直のジャンプ（この符号が正となるジャンプ）が起こり、いずれかの時点で電気機器がスイッチオフになったとき、垂直のジャンプを打ち消す符号が負のネガティブのジャンプが起こり、電気機器の使用電力は０になる、との考えによっている。なお、このことは図９に示すように実際の事例で検証済みである。

ここで、状態変数である個別電気機器の稼働状態を明確に定義しておく必要がある。ここでは、推定された個別電気機器の「機種」に属する電気機器の台数は1台であるとの単純化の仮定をおき、個別電気機器の稼働状態の状態変数を、その稼働確率と定義して、モデル化をおこなっている。しかしながら、各機種に対応する個別電気機器の台数が1台との単純化の仮定は、容易に拡張できることは言うまでもない。

以下、図３に示すフローチャートにしたがって説明すると、まずステップS1では、稼働状態推定部6が、電力計5から家屋1の総使用電力のデータを受けとる。続いて、ステップS2へ進む。

ステップS2では、時間変動算出部6aが、ステップS1で得た総使用電力に基づいて、所定時間（ここでは１分）間隔での総使用電力の時間変動１段階差の値（ジャンプ電力）を算出し、このジャンプ電力のデータを稼働電気機器機種特定部6b、対応尤度推定部6c、および稼働状態推定部6dへそれぞれ出力する。続いて、ステップS3に進む。

ステップS3では、稼働電気機器機種特定部6bが、最初から徐々に学習して記憶してきた機種（クラスタに相当）のリスト6g、各機種の規格電力（その平均値と分散）、稼働確率のデータを用いて、ステップS2で得られたジャンプ電力に基づいて、それまで記憶していた機種のものとは異なる新しい機種であるか否か、を判定する。この判定結果が、それまでにない新しい機種である、と判定（YES判定）した場合にはステップS4に進み、その判定結果が新しい機種ではない、と判定(NO判定)した場合にはステップS5に進む。
具体的には、最初の時点で、機種数が０から１種類と少ない場合には、新しい機種を加える判定基準を緩く、経過時間が長く、データ長が十分とれる時点では、情報量基準などの厳密な統計的判定基準を用いる。経過時間が短い時点では、たとえば、規格電力の大きさの10%以下のジャンプ電力であると同一の機種リスト6gのままとし、規格電力の大きさの10%より大きいジャンプ電力については、新たな機種として機種リスト6gに加える。

ステップS4では、記憶していた既存の機種（クラスタ）の数kを増やす。ステップS2で得られたジャンプ電力の絶対値を、新しい機種の規格電力（クラスタの平均に相当）として決定し、機種のリスト6gを新たに更新する。続いて、ステップS5に進む。

ステップS5では、新たに増やしたクラスタ数kに応じて機種ごとの規格電力の平均値の再計算を行う。すなわち、総使用電力の１分間隔の時系列データのうち、決められた過去のデータ長より以前のデータは削除する。具体的には、たとえば、過去１碑文1,440分あるいは１週間分10,080分以前のデータは削除し、さらに、ジャンプ電力が0でないデータを、利用するデータとして抽出する。この抽出されたデータを用いて、機種リスト6gに加えられた機種の規格電力（クラスタの平均）を再計算する。続いて、ステップS6に進む。

ステップS6では、ステップS5で再計算した規格電力の平均の計算結果にもとづき、クラスタの再分析を行って、同一の規格電力のものであるかどうか、すなわち同じ機種であるかどうかに応じて、重複する余分な機種の枝刈りを行うか否かを判定する。この判定結果がYESであれば、ステップS7に進み、判定結果がNOであれば、ステップS8に進む。

ステップS７では、既存の機種リストから重複する余分な機種が削除された場合には、余分な機種が削除された機種リスト6gの機種ごとに規格電力（平均）を再計算する。この計算結果は、機種リスト6gへ出力する。続いて、ステップS8に進む。

ステップS8では、更新された機種リストの機種ごとの規格電力の分散を求め、機種リスト6gに関するデータ（機種ごとの規格電力（平均）と分散）を対応尤度推定部6cへ出力する。続いて、ステップS9に進む。

ステップS9では、対応尤度推定部6cが、時間変動算出部6aから入力されたジャンプ電力と、稼働電気機器機種特定部6bから入力された機種リスト6gに関するデータ（機種ごとの規格電力の平均値と分散）と照らし合わせて、ジャンプ電力の大きさおよびジャンプの符号の正負からどの電気機器がスイッチオン/スイッチオフになったのか、の尤度（確率）を推定する。この尤度に関する情報は、個別電気機器稼働状態推定部6dへ出力する。続いて、ステップS10に進む。

ステップS10では、個別電気機器稼働状態推定部6dが、マルコフスイッチングモデル6eを用いて、ステップS8で得られた、現在の個別電気機器の稼働状態のもとで総使用電力の変動ジャンプから得られる個別電気機器のスイッチオン/スイッチオフの事象生起に関する尤度が入力されたとき、どのように個別電気機器の稼働状態が変化するか、を推定すし、現在の個別電気機器の稼動状態の推定を更新する。この推定結果は、カルマンフィルタ6fへ出力する。続いて、ステップS11へ進む。

ステップS11では、カルマンフィルタ6fを用いて、稼動電気機器機種特定部6bから入力された機種リスト6gに関するデータ、より具体的には、機種ごとの規格電力の平均値とステップS9でマルコフスイッチングモデル6eから入力された現在の個別電気機器の稼働状態の更新された結果から総使用電力の予測値を計算し、ステップS1で測定された総使用電力と比較し、リアルタイムで個別電気機器の機種推定、規格電力推定および稼働状態の推定を最適化する。

以上のように、個別電気機器の稼働状態推定プロセスは、大きく分けて、図３に示す領域A、B、Cで示す３つのプロセスからなる。
すなわち、同図での領域Aは、ステップS2〜S9からなり、総使用電力の時間変動１段階差データから稼働電気機器の機種を抽出して特定するプロセスに対応する。この領域Aのプロセスは、新規の個別電気機器を既存の電気機器の機種リストに新たに加えるか否かの判定し、また、既存の機種リストから同じ機種と考えられる機種を削除するために実行する。
一方、領域Bは、ステップステップS9からなり、総使用電力の時間変動１段階差、すなわちアップおよびダウンのジャンプがどの電気機器のスイッチオンとスイッチオフの事象の生起に対応しているか、の尤度を推定するプロセスに対応する。
また、領域Cは、各ステップS9〜S11からなり、測定された総使用電力の時間変動の時間間隔で刻々変化する個別電気機器の稼働状態を動的に確定するプロセスに対応する。

さらに詳細には、領域Aのプロセスは、領域A-1で示すステップS3〜S4の新しい機種を抽出するプロセスと、ステップS4の規格電力の再計算と余分な機種の枝刈りのプロセスからなる。
また、領域Cのプロセスは、ステップS10の、現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる個別電気機器のスイッチオン/スイッチオフの事象生起に関する尤度（確率）がインプットされたとき、どのように個別電気機器の稼働状態が変化するか、を、マルコフスイッチングモデル6eを用いて推定するプロセスと、ステップS11の総使用電力の時間変動の測定データから、個別電気機器の機種推定、規格電力推定および稼働状態の推定を、カルマンフィルタ6fを用いて最適化するプロセスとからなる。

図３に示す上記個別電気機器の稼働状態推定プロセスにあっては、この全体の流れは、状態空間モデルのモデリングの方法に従っている。
すなわち、上記領域A、Bで示すプロセスは、可変k -ミーンズアルゴリズムに、また領域B、Cで示すプロセスはカルマンフィルタによるモデリングの部分である。そして、両者に共通する領域Bで示すプロセスが、k -ミーンズアルゴリズムとカルマンフィルタとを接合するプロセスであって、前者の出力を後者の入力とする部分である。より詳細には、領域Aでのプロセスは、k -ミーンズモデルを、クラスタの数kを可変に拡張しEMアルゴリズム（Expectation Maximization Algorithm：期待値最大化法）と呼ばれるアルゴリズムを用いて不完全データの下で最尤推定値を求めるプロセスである。この最尤推定値は、個別電気機器の規格電力の平均と分散に関する最尤推定値である。

以上、NIALMのためのリアルタイム逆推定アルゴリズムの概要であり、上記説明では数式を使わずに概観した。
以下、これらのプロセスの詳細について、数式を用いながら説明する。

まず、領域A、Bで実行するk -ミーンズアルゴリズムについて説明する。
k -ミーンズモデルは、混合正規モデル(Gaussian Mixture Model)とも呼ばれるモデルである。サンプルが未知のk個の正規分布からの実現値であるとき個の正規分布の平均と分散を推定し、各サンプルがどの正規分布からの実現値であるのかを判定するモデルである。したがって、一般的にk -ミーンズモデルは、クラスタリングのモデルの一種である。しかし、クラスタリングに、クラスタ分けの外的基準を用いないため、教師なし学習アルゴリズムの一種として知られている。
k -ミーンズの解法には様々な考え方が取り入れられているが、本実施例では、上記のように統計的な概念がもっとも明確な混合正規モデルの立場をとり、解法としては、これも統計学の明確な理論的背景をもつ、EMアルゴリズムを採用する。

なお、以下の記号にあって、説明文章中の記号は、［数１］〜[数１１３]での記号とフォントとは、便宜上、異ならせてあるが、これらは同じものである。また、記号の上付き文字と下付き文字とが上下に重なる場合には、説明文中においては左右（下付きを最初に、その後に上付き）にずらして記載している。
混合正規モデル視点からすk -ミーンズモデルは、本実施例の場合、単純な1次元モデルである。
すなわち、m個のサンプルをx_i, (i=1,…m)と表わし、また、k個の未知の正規分布N (μ_j,σ_j ²), (j=1,…,k)と表す。ここでkは既知である。なお、説明の便宜のため、ここでは、σ_j, (j=1,…,k)は既知とし、σ_j =σ^, (j=1,…,k)で、すべて等しいとする。k=2の場合の模式図を図４に示す。

この図から明らかなように、サンプルからみると、x_iが2つの正規分布N(μ_j, σ^²), j=1,2から得られる尤度は、それぞれ、次となる。

いま、収束計算の今期ステップ目での、平均μ_j, j=1,2の推定値をμ_j ^(s), j=1,2とすると、式(1)は、sステップ目に、サンプルx_iが、2つの正規分布N(μ_j, σ^²), j=1,2から得られる尤度の推定値と考えられるので、これを同様に、l_y ^(x), (j=1,2; i=1,…,m)と表す。k -ミーンズモデルのEMアルゴリズムによる解法では、各サンプルが、未知の2つの正規分布から得られる尤度を最大化する、未知パラメータμ_j, j=1,2を、次のウェイト付平均の収束計算で求める。
式(2)の右辺は、式(1)から、μ_j ^(s), j=1,2の関数であり、μ_j, j=1,2の繰り返し収束計算となっている。このように、まず、はじめに、未知のパラメータμ_j, j=1,2の初期値μ_j ⁽⁰⁾, j=1,2を決め、式(2)にしたがって、順次、μ_j ^(s), s=0,1,2,…を計算し、収束したところで、未知パラメータμ_j, j=1,2の推定値μ_j ^, j=1,2を求める手順である。

以上の説明から、サンプルx_iが、多次元になっても、分散共分散行列の設定はあるものの、式(1)の尤度は、同じように得られるので、全く、同様に扱える。
ここで、式(1)が収束するのか、また、収束したとして、混合正規モデルの尤度を最大化しているのか、ということについて説明する。
実は、EMアルゴリズムとは、不完全データの下で、最尤推定値を与える、一般的なアルゴリズムのことである。このEMアルゴリズムは、大域最適化は保障されないが、尤度の単調(増大)性をもったアルゴリズムである。
以下では、k -ミーンズモデルを、不完全データの下での最尤推定法の問題に変換し、混合正規モデルとしてのk -ミーンズモデルを解くEMアルゴリズムを導出する。

不完全データの厳密な定義は後述する。また、不完全データを含んだモデルの尤度や最尤推定法についても、後で、一般的な解説をおこなう。ここでは、不完全データを、観測されない確率変数、潜在変数（latent variable）として捉え、まず、k -ミーンズモデルを定式化する。
混合正規モデルとしてk−ミーンズモデルのキーとなるポイントは、サンプルが、k個の正規分布N(μ_j, σ_j ²), j=1,…,kのどの正規分布からのサンプルであるかわからない、つまり、観測されていないことである。もしもこれがわかっているのであれば、通常の最尤推定法で、k個の正規分布の平均、分散の最尤推定値を求めることは容易である。通常、不完全データモデルでは、観測されないが、観測されたとしたときの確率変数を、潜在確率変数として導入してモデル化する。

そこで、1次元のk -ミーンズモデルに再びもどって、潜在確率変数を導入して、k -ミーンズモデルを定式化する。
確率変数Z_ijを、観測されないが、x_iが、第j番目の正規分布N(μ_j, σ_j ²)からのサンプルであるとき1、そうでなければ0の値をとる潜在確率変数とする。

さて、観測されない確率変数Z_ijを想定し、確率変数Z_ijが観測されたと仮定したときの完全データy_i, i=1,…,mによる表現を次で与える。

完全データ：
統計モデル：
上述したように、ここでサンプルx_iは観測されるが、Z_ij, j=1,…,kは観測されないものとしている。

通常の手順で、完全データy_i=(x_i, z_ij, j=1,…,k), i=1,…,mに対応する尤度関数を求めると以下のようになる。ここで、未知パラメータ(μ_j, σ_j ²), j=1,…,kをまとめて、仮設h=[(μ_j, σ_j ²), j=1,…,k]と表現しておく。現在の仮説を、h'=[(μ_j', σ_j'²), j=1,…,k]とするとき、この仮説のもとでの完全データの尤度は、次式 となる。
計算の簡便化のために両辺の対数をとって、全サンプルの尤度を求めると次式となる。ただし、ここで、y=(y_i, i=1,…,m)は完全データベクトルである。

この対数尤度を最大化する、(μ_j, σ_j ²), j=1,…,kの最尤推定値(μ_j^', σ_j^'²), j=1,…,kは、式(8)を(μ_j', σ_j'²), j=1,…,kに関して最大化して求められる。それぞれ以下のようになる。
平均μ_j', j=1,…,kの最尤推定値μ_j^', j=1,…,kは、式(8)をμ_j', j=1,…,kに関して微分し、1階の条件から、μ_j^', j=1,…,kを求めれば、次式が得られる。

また、分散σ_j'² , j=1,…,kの最尤推定値σ_j^'² , j=1,…,kは、同様に、式(8)を、σ_j'² , j=1,…,kで微分し、1階の条件を求めると、次式が得られる。
これを解くと、次式が得られる。
したがって、最尤推定値(μ_j^', σ_j^'²), j=1,…,kを式(8)に代入すれば、最大対数尤度l^*を求めることができる。

ここで、
を代入すると、式(12)は以下のように表せる。

さて、本実施例では、クラスタの数kが、固定されているのではなく、可変の場合を考える。最適なクラスタ数k*の決定には、赤池の情報量基準(AIC)を利用することができる。
このAICを用いて、最適なクラスタ数を、AICを最小にするk*とすることができる。AICは最大対数尤度をl*、パラメータ数をpで表すと、次式で定義される。

k -ミーンズモデルでは、クラスタの数kを１つ増やすと、平均と分散の2つのパラメータが増加することに注意する。また、クラスタの数kに依存して、最尤推定値(μ_j^', σ_j^'²), j=1,…,kも最大対数尤度i*も、また、AICも変わるので、これらをAIC(k)のように、kの関数として表すことにしよう。
最適なクラスタ数k*を求めるアルゴリズムの表現としては、順次、kを1から増やしていき、AIC(k+1)>AIC(k)となるkを、最適なクラスタ数k*とすればよいことになる。

具体的に、AIC(k+1), AIC(k)を求めると、次式(15)、(16)が得られる。
したがって、次式が成り立つkを求めることになる。

不完全データ最尤推定法では、EMアルゴリズムやこれを点対集合写像に拡張したGEM(Generalized Expectation Maximization)アルゴリズムを、不完全データの下での重力モデルの推定やIPFP(Iterative Proportional Fitting Procedure)の拡張に適用する。まず、不完全データの定義、GEMアルゴリズムの特徴づけについて説明する。
確率変数YとXを考える。データとは、これら確率変数の実現値であり、小文字y, xで表す。一般に、確率変数yとxとの間に、多対一の写像hによる対応関係があるとき、yを完全データ、xを不完全データと定義する。すなわち、
データで表現すると、同じことであるが、完全データyと不完全データxとの間に以下の関係を想定したことになる。

この定義は、一般的、抽象的な定義であるが、理解を確実にするために、潜在確率変数Z_ij, i=1,…,kが、不完全データの例となっていることを、確認しておく。
サンプルとなるデータベクトルをx=(x_i; i=1,…,m)、また、潜在確率変数に対応する潜在変数データベクトルをzと表そう。zはその定義より、z=(z_ij; i=1,…,m, j=1,…,k)である。
完全データyは、xとzの組、
である。ここで、多対一の写像hを次式のように定義する。
実際、写像hが多対一となっていることは、zのとるどんな値に対しても、zと対をなしているxの値が定まれば、写像hの像は常にxとなっていることからわかる。これは、(x, z)の2次元平面をx軸に射影したものを想像するとわかりやすい。

さて、観測されるのが、確率変数xのみであるとき、別言すれば、観測されるのが、不完全データxのみであるとき、最尤推定法では、この不完全データを生成する尤度を最大にする未知パラメータの値を求めて、当該パラメータの最尤推定値とする必要がある。
その方法は、あるモデル（パラメータ）の下で、不完全xデータを生成する確率（尤度）を、完全データyのモデルに逆もどって、計算することである。
つまり、写像の像xに対応する原像が、像xを生み出す完全データyの集合であることに注意すると、不完全データxを生み出す確率(尤度)は、原像h^-1(x)に含まれる完全データyの起こる確率である。別言すれば、原像h^-1(x)に含まれる、個々の完全データyが生起する確率（尤度）の和が、不完全データxが起こる尤度となる。以上を数式で表現すると、以下のようになる。

まず、不完全データxを生起させる完全データyの集合を次で定義する。
不完全データxが起こる確率(尤度)をg(x|Φ)で表せば、以下のようになる。
したがって、不完全データの最尤推定法は、式(23)の尤度を最大にするパラメータΦを求めることである。

一般に、不完全データの尤度g(x|Φ)をΦに関して最大化する問題の解法は簡単ではない。この解法としては、EM (Expectation Maximization)やGEM (Generalized Expectation Maximization)といわれる方法が考案されている。
これらの方法の概要は以下の通りである。
まず、不完全データの対数尤度をL(Φ)で表す。
次に、不完全データxおよびパラメータΦが与えられたときの完全データyの条件付き確率をk(y|x, Φ)で表す。

また、Q(Φ'|Φ)とH(Φ'|Φ)を、それぞれ次のように定義する。

以上の定義から、対数尤度L(Φ')が、次式で表される。
なぜなら、一般的に、E[h(x)|x]=h(x)が成立することで、
が成り立つからである。

EMアルゴリズムは、次の2つのステップからなる、第ステップにおけるΦの推定値Φ⁽ⁿ⁾から、第ステップの推定値Φ⁽ⁿ⁺¹⁾を計算する、Φ⁽ⁿ⁾→Φ⁽ⁿ⁺¹⁾, n=0,1,2,…の繰り返し計算法である。ここで、新しいステップであるか否かを判定する期待値ステップ（E-ステップ）と、平均を評価する最大化ステップ（M-ステップ）とでは、それぞれ以下の計算を行う。
EMアルゴリズムの定義：
E−ステップ:
M−ステップ:

一方、GEMアルゴリズムは、EMアルゴリズムのMステップを、次の点対集合写像Mに置き換えたものである。

GEMのM−ステップ:
ここで、点対集合写像Mは、次のように定義される。

なお、次式が成立することが知られており、
これを用いて、GEMの単調(増加)性が証明されている。
したがって、対数尤度L(Φ)が、上に有界であれば、L(Φ⁽ⁿ⁾)は、どこかに収束する。ただし、関数値L(Φ⁽ⁿ⁾)の収束が、ただちにΦ⁽ⁿ⁾の収束を保証するものではない。

以上の準備のもとで、k-ミーンズモデルのEMアルゴリズムによる解法を構成する。
まず、EMアルゴリズムでは、E−ステップの計算が必要である。k -ミーンズモデルのパラメータはh[(μ_j,σ_j ²), j=1,…,k]であるので、Φ, Φ'の代わりに、h, h'の表記を用いる。
また、E−ステップの条件付き期待値は、条件付き対数尤度の条件付き期待値となっており、複雑な計算と捉えられるかも知れない。しかし、条件付き対数尤度の条件となっているパラメータは、条件付き期待値の計算に何の影響も与えないことに注意すると、式(26)の計算は、と不完全データxの下での完全データyの期待値を計算するにすぎないことになる。

それでは、k-ミーンズモデルに対するQ(h'|h)、すなわち、式(26)の対数尤度の条件付き期待値Q(Φ'|Φ)を求めよう。Q(h'|h)は、式(8)を参照して、次のように表すことができる。
GEMアルゴリズムの特徴は、Q(h'|h) ? Q(h|h)なる任意のhを選んでも、尤度がL(h') ?L(h)の単調性をもつことである。

前述のように、EMアルゴリズムとは、Eステップで式 (36) を計算し、Mステップで、h⁽ⁿ⁺¹⁾=argmax_hQ(h|h⁽ⁿ⁾)として、{h⁽ⁿ⁾, n=1,2,…}の系列を生成し、収束点を不完全データの最尤推定値とする方法である。
さて、Q(h'|h)の式(36)をみればわかるように、Eステップにおいては、次の計算が問題となる。

ここで、式(37)の計算方式について、説明する。
式(37)は、事後確率の形をしている。すなわち、不完全データx_iを観測したもとで、x_iがどの正規分布からのサンプルかを判別するインジケータ確率変数Z_gの事後確率を求めている。この事後確率はわからないものの、インジケータ確率変数Z_gが与えられたもとで、不完全データx_iが起こる尤度は明確に計算できる。
つまり、ベイズの定理を用いることで、式(37)の事後確率を、事前分布と尤度の積に比例するものとして求める。

そこで、Z_g=1となる事象をΩ_yと表すと、次式が成立する。
一方、すべてのj≠j'に対して、Ω_yとΩ_{y '}は、明らかに、排反であるので、
また、
が成立する。したがって、期待値と確率との間に次の関係が成り立つ。
これを用いると、ベイズの定理より、
が得られる。

事前確率P(Z_y)を、不完全データx_iが、時間順序tにしたがって、観測されたと仮定し、t期直前のX_t-1,…,X₁までの観測に依存する事前分布に拡張しよう。
πj|i-1=P(Z_y|X_i-1,…,X₁)とおくと、式(43)は、次式で表される。
とおいて、本実施例で用いるk-ミーンズの解法アルゴリズムとしてのEMアルゴリズムは、次のようにまとめることができる。

Eステップ：
Mステップ：

式(46)および式(47)から、このEMアルゴリズムは、h=[μ, σ²]を入力として、h'=[μ^', σ^'²]を生成する、繰り返し計算法となっていることがわかる。
以上の説明では、複雑さを避けるために、クラスタの数kが固定されたケースを考えてきた。しかし、本実施例で構築したNIALMのためのリアルタイム逆推定アルゴリズムでは、クラスタの数kが可変のk -ミーンズモデルを採用している。

次に、クラスタの数kが可変のk -ミーンズモデルについて説明する。
本実施例では、観測された1分おきの総使用電力にもとづいて、逐次的に電気機器の機種を同定し、個別電気機器の機種の数が変動するモデルを目標としている。総使用電力の1分時間階差をベースに、どのような機種の電気機器があるのか、また、それらの電気機器の稼働状態はどうなのか、をリアルタイムに推定するアルゴリズムが目標である。
本実施例では、その電気機器の機種の同定にk -ミーンズモデルを適用する。通常のk -ミーンズモデルでは、クラスタ数kは、事前に与えられた数に固定されている。しかし、電気機器の種類は、k -ミーンズモデルによるクラスタに対応するから、当然、k -ミーンズモデルのクラスタ数kも可変となる必要がある。

応用面からみて、なぜ、逆推定のアルゴリズムが、個別電気機器の機種の同定まで行う必要があるのか、について以下に説明する。
実用化の場面で、日本国内で3000万世帯にもおよぶ個々の家庭にHEMS(Home Energy Management System)を導入する場面を想定してみる。逆推定のアルゴリズムをこれらの家庭に実装する場合、それぞれの家庭が持っている電気機器を種類まで特定し、事前情報として持っておかないと、アルゴリズムが動かないので、大変な実装の手間とコスト、時間がかかる。より実用化の段階に近ければ近いほど事前情報にたよらず、アルゴリズム自体が学習し、初期設定にたよらなくてよいシステムが望まれる。本実施例もそのような実用化のレベルに近いレベルでのシステム化を想定している。

すでに、式(15)〜(17)において、AICを用いた最適クラスタ数の決定方式について記述した。
理論的には、これで十分のように思われるが、実用化を想定したモデリングでは、これでは不十分である。
なぜなら、AICの定義に中に、大標本理論での極限の概念があり、一体、どの程度のサンプル数であれば、不偏推定量として、どの程度の誤差の範囲に収まるのかの、議論がないからである。
実際、アルゴリズムがスタートして直後は、データ数（サンプル数）が少なく、アルゴリズムも学習が進んでおらず、個別電気機器の機種の付け加えや、削除などをAICの基準で行おうとすると、すぐに破綻する。この問題は、一般に、コールドスタートの問題とも呼ばれ、リコメンデーションシステムなどのシステム設計でも、新規のユーザーについては、利用履歴が少なく、システムがうまく作動しないことが起こる問題として知られている。

そこで、本実施例では、コールドスタート問題を考慮し、クラスタ数の決定については、アドホックな、実用的な方法を採択した。
ここでは、クラスタ数は、個別電気機器の機種の数を表し、また、クラスタを構成する正規分布の平均が、個別電気機器の規格電力に対応している。以下に、本実施例で採択した可変クラスタ数の決定方法を示す。
本実施例ではまず、総使用電力の1分時間間隔階差の電力ジャンプ量に着目する。

図３の領域A-1での、ジャンプ量から新たな電気機器の機種を同定し、既存の機種リストに加えるプロセスは、
1) ジャンプ量が0に近い場合
個別電気機器の稼働状態に何も変化なしとして、次の期に移る。
一方、
2) ジャンプが0ではない場合、
a) ジャンプ量が既存の電気機器の規格電力に近い場合は、既存の電気機器がオンまたはオフしたと判定し、
b) ジャンプが既存の電気機器の規格電力から離れている場合は、新たな電気機器の機種が特定されたと判定し、新たな電気機器の機種の規格電力をそのジャンプ量に設定してこのジャンプ量をともなった新たな機種を既存の機種リストに加える。

このように、上記の新しい電気機種の導入判定プロセスでは、1) 既存の電気機器の機種リストに新たな機種を加えるか、あるいは、2) 何も加えないか、のいずれかになる。すなわち、新しい機種と判定する基準としては、ジャンプ量が、規格電力の分散のある一定比率を超えた場合に、新しい機種として判定し、また、そうでない場合を、既存の機種の電力使用が規格電力の通常の変動範囲内に入っていると判定して、既存機種リストに変更を加えない。

領域A-1のプロセスで決定された電気機器の機種リストをもとに、EMアルゴリズムを実行し、個別電気機器の規格電力の更新、再計算を行う。
これが、次の電気機種の重複削除（枝刈り）のプロセスである。
すなわち、領域A-2のプロセスにおいて、 規格電力およびその分散の更新と機種の重複削除を行う。

1) 領域A-1のプロセスからの得られる電気機器の機種リストとその規格電力および規格電力の分散をもとにEMアルゴリズムを1ステップ計算することで、個別電気機器の規格電力および規格電力の分散の再計算による更新を行う。
2) 規格電力の再計算の結果、規格電力が同一とみなされる機種を1つの機種に統合する。すなわち、重複機種の削除を行う。
このように、領域A-2のプロセスは、EMアルゴリズムによる平均、分散の再計算、更新をおこなうプロセスである。ここでも同様に、同一の機種に属すかどうかの判定は、規格電力の分散に対する大きさを判定基準としている。

さて、EMアルゴリズムは、本来、繰り返し収束計算を何回もおこなって、収束した値を、平均、分散の推定値とするアルゴリズムである。しかし、1分間隔で入力されてくる総使用電力やそのジャンプ量を用いて、1分間隔ごとに、各期で収束するまで、収束計算を繰り返すとすれば、計算コストが非常に高くなる。
そこで、本実施例では、尤度の単調性の特徴をもつEMアルゴリズムの特性を利用し、1回のみ、EMアルゴリズムの繰り返し計算を行なって、次のプロセスに引き継ぐことにしている。
同一のこれまでに登録された規格電力の値に近ければ、その規格電力に対応した電気機器がオン、オフしたと考える。既存の規格電力の値に近いか否かの判定は、規格電力に対応した分散の大きさで判定する。
なお、上記Mステップにあっては、用いる過去のデータの長さ、逆戻る長さは、前記過去のデータより少ない長さに指定する。
以上で、クラスタ数が可変の可変k-ミーンズ アルゴリズムの説明を終えたので、次に、NIALMのための逆推定アルゴリズムの動的表現について説明する。

NIALMのためのリアルタイム逆推定アルゴリズムの目的は、実際の家庭環境で、総使用電力の1分時間間隔の変動データ（ジャンプ量やそのアップまたはダウン方向）から、リアルタイムに、どのような電気機器が存在し、どのような稼働状態にあるかを推定することであった。そのためには、時系列で変動するデータをリアルタイムに動的に処理するシステムが必要である。さらには、時系列データのリアルタイムの処理のみではなく、リアルタイムの制御にまで結びつける必要がある。

現在のところ、時系列データをリアルタイムに動的に処理しながら、制御にまで結びつける方法として、もっとも有力な方法として、状態空間モデリングと呼ばれている方法が知られている。
この状態空間モデルを用いる方法は、カルマンフィルタとも呼ばれる。その基本的考え方は、潜在変数である状態変数の時間的遷移を表現する状態方程式と、観測されない状態変数を観測データに結びつける観測方程式とから構成される。
本実施例では、状態方程式が、総使用電力の1階階差のジャンプからどの電気機器がオン/オフしたのかといった統計的推定にかかわるモデルの記述式としている。

状態空間モデリングを考える際の重要な事項は、1) 状態変数をどのように定義するか、2) 状態がどのように時間的に遷移するかを記述する状態方程式、そして、3) 観測されない潜在変数としての状態変数を、観測変数に結びつける観測方程式、である。
さらに、4) 制御と関連づける場合は、状態変数以外に、制御者や観測者が、状態方程式で表される動的システムに、政策変数、あるいは、制御変数として、入力を加える必要がある。
本実施例のNIALMのためのリアルタイム逆推定アルゴリズムは、この最後の入力変数を含んだモデル化を行ったことを特徴とする。すなわち、本実施例での入力変数を含んだモデル化では、旧来の制御変数の意味ではなく、動的システムの状態を変える外的なショックと捉えたモデル化を用いる。もっと広義に一般化すれば、外部からの学習入力なども、同じような扱いが可能になる方法である。

ここで、カルマンフィルタの表現の標準形を示す。
まず、t期における状態変数ベクトルをx_i、入力ベクトルをu_t、観測ベクトルをy_tとおく。状態方程式の誤差項の確率変数ベクトルをW_tとおき、観測方程式の誤差項の確率変数ベクトルをv_tと表す。
標準形では上述の状態方程式、観測方程式は、それぞれ以下のように表す。
ただし、F_t, G_t, H_tは、既知の係数行列である。

ここで、状態方程式の誤差項、また、観測方程式の誤差項は、異時点間で独立である。また、状態方程式の誤差項と観測方程式の誤差項は、互いに独立とする。
以上の仮定は、次で表現される。
ただし、δ_n'は、t=t'のとき1、それ以外のとき0のクロネッカーのデルタである。また、式(52)のOは、適当なサイズのゼロ行列である。また、Q_t, R_tは既知の共分散行列である。

ここで、入力ベクトルの時間添字について、後の記述にも関連するので、注意しておく。入力ベクトルの時間添字がt-1となっているが、式(52)は、t-1時点の状態変数を与件とし、入力が加わり、t時点の状態変数を生み出すことを示しているに過ぎない。つまり、入力の時点sは、t-1 < s < tであればよく、t-1となっているのは、便宜的な表記であって、tでもよい。
また、誤差項の時間添字がtとなっているが、これもt-1時点の状態変数に誤差項が加わって、t時点の状態変数を生み出すことを示すためである。入力ベクトルは、観測誤差を含まないと想定しているので、入力の時点sは、誤差の加えられる時点の前でも後でもよい。このことから、後の記述では、入力の時点をtとしているが、本質的に式(52)と変わらないことを付記しておく。

ここで、次の条件付き確率の表記法を導入する。そうすると、カルマンフィルタの再帰的改定の表現を簡明化することができる。
すなわち、上記各式にあって、右肩の添字は、1期前までの観測値が与えられたときの条件付き期待値を表し、右肩の添字は、今期までの観測値が与えられたときの条件付き期待値を表している。

同様に、状態変数ベクトルの共分散行列の条件付き期待値にも同様の表記を用いる。
これを、略記して、次のように表す。
P_t ^-,P_t-⁺ は、いずれも状態変数ベクトルの誤差共分散である。

初期ステップは次のように与える。
ここで、入力の系列、u₀,u₁,…u_t,…が問題であるが、入力の系列は、誤差なく外部から与えられていると仮定している。

カルマンフィルタの核である再帰的(Recursion)な更新プロセスは、次のように表現される。
式(63)は、カルマンゲインと呼ばれる項である。

以上の方程式体系にしたがったカルマンフィルタの計算手順は次のようになる。
1) まず、式(59)および(60)にもとづいて、初期値(x₀^⁺, P₀ ⁺)を決める。
2) 次いで、式(61)および式(62)にもとづいて、状態変数ベクトルの予測値x_t^^-と、誤差共分散P_t ^-を１期前までの観測値にもとづいて求める。
3) 式(63)のカルマンゲインを計算する。
4) 式(64)により、状態変数ベクトルの予測値を今期の観測値にもとづいて改定する。
5) 式(65)により、状態変数ベクトルの誤差共分散を改定する。
以上が、状態空間モデリングの考え方、カルマンフィルタの標準形、そして、カルマンフィルタの再帰的計算の手順についての説明である。

次に、上記カルマンフィルタを実際にNIALMのためのリアルタイム逆推定アルゴリズムの構成に適用する。
NIALMのためのリアルタイム逆推定アルゴリズムは、総使用電力の1分間隔の時間変動データから、個別電気機器の機種（それぞれ規格電力を有する）の同定・推定と、各電気機器の稼働状態の動的変動をリアルタイムに把握することをねらったアルゴリズムである。本実施例では、状態空間モデリングの方法を個別電気機器の稼働状態の動的推定に適用するとともに、これを前述した個別電気機器の規格電力の推定アルゴリズムと統合することで、逆推定アルゴリズムの目的を実現する。
その最も特徴的な点は、個別電気機器の規格電力推定アルゴリズムによる推定結果を、カルマンフィルタの入力として活用し、個別電気機器の稼働状態の状態遷移を表す状態方程式に組み込んだ点である。

具体的な個別電気機器の稼働状態の動的推定を扱うため、以下の記号を導入する。
Pw_t: 第t期の総使用電力
Pw_t ^k: 第k種電気機器の第t期の電力使用
P^k: 第k種電気機器の規格電力
(σ²)^k: 第k種電気機器の規格電力の分散
P_t^^k: 第k種電気機器の第t期での規格電力の推定値
(σ^²)_t ^k: 第k種電気機器の第t期での規格電力の分散の推定値

観測値は、総使用電力の1次元の1分間隔観測データである。

y_t=Pw_t, t=1,2,…

なお、カルマンフィルタの標準形に記法を合わせるために、y_tの記号を用いている。

さて、状態変数であるが、NIALMのためのリアルタイム逆推定アルゴリズムの目的が、個別電気機器の稼働状態の動的推定であるから、当然、状態変数ベクトルは、各電気機器の稼働状態を表したものとする必要がある。これをより明確にするために、稼働状態がオン/オフをとる潜在確率変数を以下のように定義する。
第種の電気機器の稼働状態を表す確率変数をX_t ^kとする。状態確率変数X_t ^kは、以下のように定義される。
ここで、状態確率変数の期待値X_t ^kをとると、
となり、期待値x_t ^kは、第種の電気機器の稼働確率となることが分かる。そこで、あえて表記の誤用を行って、以下では、状態確率変数の期待値x_t ^kを状態変数とみなして、状態変数ベクトルは、各電気機器の稼働確率を指しているものとみなしていく。
ただし、例外的に、機種0は、総使用電力を表し、状態変数に加えておく。

以上、まとめると、以下となる。
x_t ⁰: 総使用電力（観測値に同じ）の期待値 X_t ⁰=Pw_t
x_t ^k, (k=1,2,…,k_t): 第種電気機器の稼働確率
ただし、k_tは、t期までに推定された個別電気機器の種類の数
Ω_t: t期までに推定された全電気機器の機種のリスト
すなわち、Ω_t={1,2,…,k_t}
B_t: t期までに推定された全電気機器の規格電力P^k、分散(σ²)_t ^kの推定値のリスト、すなわち、
したがって、K_t=|Ω_t|である。ただし、|S|は、集合Sの要素の数を表す。

本実施例の逆推定アルゴリズムの最大の特徴は、総使用電力の1分間隔の時間階差であるジャンプを、カルマンフィルタの入力変数に取り入れている点である。
ただし、ジャンプをそのまま入力にしているわけでではなく、ジャンプが既存の保有電気機器の集合Ω_tの中のどの機器の稼働状態の変化から生じたのかを判定し、各状態変数への入力ベクトルに、k-ミーンズモデルを用いて変換している点が、カルマンフィルタの標準形と異なるところである。
まず、総使用電力の1階階差としてのジャンプを、スカラーの入力変数とし、総使用電力の状態変数X_t ⁰に対応するスカラーの入力変数とする。
通常、入力変数には、誤差がないものとして扱うので、u_t ⁰=E(u_t ⁰)であり、次式が成り立つ。

以後、入力変数については、期待値と区別せずに扱っていく。
次いで、k_t種類の電気機器に対応する状態変数への入力変数u_t ^k, k=1,…,k_tを定義する。
これらは、u_t ⁰のジャンプが、1,…,k_tのどの機器の稼働状態の変化から発生したのかを表すk_t個の電気機器の種類上の確率分布として定義される。これらの入力変数u_t ^k, k=1,…,k_tについては、この後、形式的な定義を述べるが、その意味については、改めて詳述する。当面は、k_t個の状態に対応する入力ベクトルであるとしておく。
これらk_t+1個の入力変数をまとめてベクトルu_tと表すことにする。

以上の準備のもとで、状態方程式を構成する。そのためには、x_t ^k, (k=0,1,…,k_t)のk_t+1個の状態変数について、x_t ^kが、x_t-1 ^kとu_t ^kから、どのように生成されるかを記述していけばよい。
1) u_t ⁰およびx_t ⁰の決定
この定義より、状態変数x_t ⁰の系列は、総使用電力Pw_t（観測値y_t）の系列と同一になる。式(71)から式(72)は、総使用電力のジャンプに着目してモデル化することを具体的に表したものである。

2) u_t ^k, (k=1,…,k_t)およびx_t ^k, (k=1,…,k_t)の決定
(a)nzu_t ^k, (k=1,…,k_t)を次で決める。
ただし、
Φ(x|μ,σ):平均μ,標準偏差σの正規分布N(μ,σ²)の密度関数
|u_t ⁰|:u_t ⁰の絶対値
π_t ^k, k=1,…,k_t：k_t個の仮説N(μ_t ^k,(σ²)_t ^k)に対する事前分布 (Cf.式(43))
式(73)は、総使用電力のジャンプがどの機種から発生したかを判別する確率を表している。

(b) 状態変数 x_t ^k, k=1,…,k_tの決定（状態方程式のうちk=0を除く、k=1,…,k_tの方程式）
次式は、状態方程式のうち、k=0を除く、k=1,…,k_tの方程式を表している。
ここで、u_t ^k, (k=1,…,k_t)を次式で定義する。
次式は、カルマンフィルタの標準形に帰着する。

3) 観測方程式の決定
P_t ^kx_t ^k, (k=1,…,k_t)は、第k種の電気機器の規格電力P_t ^kと稼働確率x_t ^kとの積であるから、第k種の電気機器のt期の電力使用の予測値となっている。また、Σ_k=1 ^ktP_t ^kx_t ^kは、総使用電力の予測値となる。これは、必ずしも、総使用電力の観測値y_tと一致しない。v_tは誤差項である。
以上のように、式(71)から式(72)が状態方程式を構成し、式(76)が観測方程式となっているので、状態空間モデリングによるカルマンフィルタの標準形に、NIALMのためのリアルタイム逆推定アルゴリズムを帰着できる。

より正確にいえば、カルマンフィルタの標準形は、再掲すると、次の線形の形をした標準形であった。
本項での記述から明らかなように、総使用電力u_t ⁰が、直接、上の式(48)に入るのではなく、式(73)の非線形の変換を経た、k_t個の入力u_t ^k, (k=1,…,k_t)が、式(48)で表現されるk_t本の状態推移方程式に代入されるので、本質的には、非線形のカルマンフィルタとなっている。

式(48)および式(76)は、k -ミーンズモデルとカルマンフィルタとの統合の核となる部分である。式(73)は、図３中の領域BAの部分に相当する部分で、k -ミーンズアルゴリズムからの出力を、カルマンフィルタの入力にどのように変換するか、を計算するプロセスである。
一方、式(74)は、図３の領域C-1のマルコフスイッチングに相当する部分であって、上記入力によって、状態変数をどのように変化させるか、を計算するプロセスである。
次に、これらの計算プロセスについて詳しく説明する。

まず、k -ミーンズの出力をカルマンフィルタの入力に変換する式(73)から説明する。
期でのk-ミーンズモデルの計算フローは、電気機器の新しい機種を既存リストに加えるかどうかを判定し、次いで、その新しい機種リストのもとで、EMアルゴリズムによって、規格電力および規格電力の分散を再計算し、重複する機種があれば削除する、といった流れである。
したがって、図３の領域Aのブロックからの出力は、t期での最新の個別電気機器の機種数、個別電気機器の規格電力および規格電力の分散の推定値である。k -ミーンズモデルの言葉でいえば、t期での最新の機種数k_t個の正規分布の平均と分散である。NIALMの用語では、個別電気機器の規格電力P^kの推定値P_t ^kと、その分散(σ²)^kの推定値(σ²)_t ^kである。
そこで、総使用電力のジャンプu_t ⁰が、k_t個の正規分布N(μ_t ^k, (σ²)_t ^k), k=1,…,k_tの中の、どの正規分布からのサンプルか、の判断をおこなうことができる。

言い換えると、総使用電力のジャンプu^t0が与えられたとき、そのジャンプu_t ⁰が、k_t個の電気機器のなかの、どの電気機器の稼働状態の変化から生じたのかを計量的に判断することができるということである。
総使用電力のジャンプu_t ⁰は、正負があるが、規格電力の大きさは正で分散も正であるから、u_t ⁰の絶対値|u_t ⁰|をとり、k -ミーンズモデルの考え方を用いて、式(44)にしたがって、求めたものである。式(73)を再掲するが、正規化していることを明記するため、式(78)に示すように、u_t ^kを、nzu_t ^k, k=1,…,k_tと表記している。
ただし、
Φ(x|μ,σ):平均μ,標準偏差σの正規分布N(μ,σ²)の密度関数
|u_t ⁰|:u_t ⁰の絶対値
π_t ^k, k=1,…,k_t：k_t個の仮説N(μ_t ^k,(σ²)_t ^k)に対する事前分布 (式(44)を参照)

nzu_t ^k, (k=1,…,k_t)は、事前分布π_t ^k, (k=1,…,k_t)のとき、u_t ⁰が観測された場合の事後分布となっている。

さて、t期の総使用電力のジャンプu_t ⁰が観測されたとき、既存の個別電気機器の稼働状態をどのように変えていくべきか、を考えよう。
本実施例での逆推定アルゴリズムの開発では、短い時間間隔で、総使用電力を測定すれば、2つ以上の電気機器が同時にオンしたり、同時にオフしたりすることが少なくなるはずだ、との前提に立っている。また、同様に、複数の電気機器が、同時にオンとオフになることも少なくなると想定している。
このような想定のもとでは、期の総使用電力のジャンプu_t ⁰は、k_t種類の機種のいずれかでオン/オフの変化が起こったと考えるのが自然である。
したがって、次のようになる。

問題は、その変化が、k_t個の電気機器の、いずれの電気機器の変化か、を推定することである。その推定のために、前述したように、k-ミーンズモデルの考え方を適用している。その結果が、nzu_t ^k, k=1,…,k_tであり、u_t ⁰が与えられたとき、それがどの電気機器のオン/オフの変化から生じたかの事後確率であった。
さて、u_t ⁰が与えられたとき、k_t個の電気機器のどの電気機器に変化が起こったのかの確率が、nzu_t ^k, k=1,…,k_tで推定されたとしても、それが状態変数にどのような変化をもたらすかは、必ずしも明らかではない。
それは、変化が起こる前のt-1期の状態変数の値x_t-1 ^kによって、状態変数の変化の仕方が異なるからである。t-1期の状態の違いによって、t期の状態遷移の仕方が異なるモデルのことを、一般に、マルコフスイッチング(Markov-Switching)モデルと呼ばれている。
本実施例におけるNIALMのためのリアルタイム逆推定アルゴリズムでは、状態遷移に、このマルコフスイッチングモデルの考え方を適用している。

以下に、これについて説明する。
まず、u_t ⁰>0の場合からみていく。図５に、u_t ⁰>0の場合の状態変化をまとめた表を示す。
u_t ⁰>0あるから、k=1,…,k_tのいずれかの機器にオンの変化があったと考えられる。しかし、この変化は、t-1期の状態変数x_t-1 ^kの値によって異なる。

その変化は、一般的に、図５の表に示すように、推移確率の形に表現できる。すなわち、t-1期に、x_t-1 ^kの値が0の場合に、t期に、x_t ^kが0,1に変化する推移確率を、それぞれ、P₀₀,P₀₁で示している。同様に、t-1期に、x_t-1 ^kの値が1の場合、t期に、x_t ^kが0,1に変化する推移確率は、それぞれ、P₁₀,P₁₁で示してある。
さて、u_t ⁰>0の場合、いずれかの電気機器がオンに変化したと考えられるから、すでに、t-1期で、x_t-1 ^kの値が1の機種、すなわち、すでにオンの機器では、変化が起こらなかったと考えるべきで、そのことを表現するために、図５の表では、P₁₁=1,P₁₀=0と記入してある。

一方、t-1期に、x_t-1 ^kの値が0の機種では、オンの変化が起こる可能性、すなわち、t期に、x_t ^kが1に変化する可能性があるので、この確率P₀₁を、nzu_t ^kで、評価しようというのが、本モデルの発想である。
この発想にしたがって、推移確率を、次のように想定する。
とすると、マルコフチェインより、x_t ^kが状態1をとる確率、すなわち、t期で、機種kの電気機器がオンになっている確率は、次式となる。

いま、確率変数の期待値は、確率を表すことに着目すると、次式(82)〜(85)が成り立つ。

結局、以上の関係を用いると、式(82)は、状態方程式となって、
となるので、改めて、入力変数を
とおくと、次のように線形標準形として表現できる。

同様に、u_t ^k<0のとき、以下となる。図６にu_t ⁰<0の場合の状態変化をまとめた表を示す。
図６の表の記載したように、推移確率を以下のように想定する。
マルコフチェインより、x_t ^kが状態0をとる確率、すなわち、t期で、機種kの電気機器がオフになっている確率は、次となる。

一方、状態方程式は
となるので、改めて、入力変数u_t ^kを
とおくと、次のように表現できる。

一方、u_t ⁰=0のときは、
すなわち、推移確率はP₀₀=1, P₁₁=1となる。
したがって、状態方程式は、次式である。

以上で、カルマンフィルタの状態方程式を導くことができる。マルコフスイッチングモデルでは、状態方程式が、単一の方程式の形式で表現きるのではなく、期の状態変数の状態に応ずると同時に、総使用電力のジャンプが、正、負、または、0の場合に応じて、複数の方程式の形式で表現されている。
さらに付け加えれば、通常、カルマンフィルタでは、状態変数ベクトルの次元は、変動せず、固定されている。しかし、本実施例の逆推定アルゴリズムでは、個別電気機器の機種と規格電力自体を、時間が経過する途上でリアルタイムに推定していくことになり、機種数が時間の経過とともに変化する。したがって、状態方程式を記述する状態変数ベクトルの次元が時間推移によって、変わることになり、状態変数ベクトルの次元が変化する際の、変数の引き渡しなどの未開拓の課題も出てくる。

次に、マルコフスイッチングモデルで記述された状態方程式と総使用電力の観測値を用いて、カルマンフィルタによる状態変数の最適推定の再帰的更新プロセスを説明する。
このプロセスは、図３の領域C-２に相当する部分である。本モデルの状態変数は、推定された個別電気機器の稼働確率であるから、総使用電力の観測値を得て、随時、個別電気機器の稼働状態の最適推定を更新していくことになる。
状態変数の再帰的な最適更新アルゴリズムは、カルマンフィルタの核である。その計算プロセスは、カルマンフィルタの標準形について、式(95)から式(65)としてすでに記述されている。

カルマンフィルタは、最小二乗法の2つの意味での拡張になっている。
第1は、通常の最小二乗法が、未知パラ―メータとしているβを、確率変数として扱い、その最少分散推定量を求めている点であり、初期時点で事前にその共分散に関する事前情報をもっている点である。
第2は、最小分散推定法において、逐次、新たな観測値が入ってきたとき、それまでのβの推定値をどのように更新していくかの論点であり、新たな観測値がもつ、それまでの推定値が張る空間への直交成分に対する寄与部分を既存の推定値に加えて、逐次、更新を行っていく点である。この第２の論点が、カルマンフィルタでは、カルマンゲインと呼ばれる項を用いた状態変数の更新アルゴリズムに表現されている。
では、順次、その計算プロセスを行列形式に変形してみていく。

式の誤差共分散R_t, t=0,1,2,…が外生的に与えられる。状態変数とその共分散行列x₀^⁺,P₀ ⁺は、t-1以後、順次、内生的に改定されていくのに対し、Q_t, R_tは、すべての期で外生的に与えられる。
次に、状態変数の右肩についたとの添字の区別が重要である。それらを再掲すると、
であり(式(53)および(54)を参照)、x_t^-は、t-1期までの情報をもとに、t期の状態変数x_tを推定した推定値である。一方、x_t^⁺は、期までの情報を用いて、t期の状態変数x_tを推定した推定値である。
カルマンフィルタの期のアルゴリズムは、x_t-1^⁺, P_t-1 ⁺が与えられたもとで、入力u_t-1を得て、順次、計算をおこない、x_t^⁺, P_t ⁺を出力し、次のt+1期の計算につなげていく、逐次計算プロセスである。まず、入力u_t-1をもとに、状態方程式によって、x_tの予測値x_t^-を次で求める。（式（48)を参照）

本逆推定アルゴリズムを具体的に行列形式で書き下すと、F_t-1, G_t-1 は、単位行列であるから、次のようになる。
次に、状態方程式の予測式(62)にしたがって、t-1期までの情報にもとづいた、状態変数x_tの共分散P_tの推定値P_t ^-を計算する。
さらに、次式によりカルマンゲインを計算する。
ただし、H_tは、t時点までに推定された個別電気機器の規格電力からなるベクトルであり、具体的には、以下のように表現される。

次に、t期の観測値、すなわち、総使用電力Pw_tを観測し、t-1期までの情報を用いたt期での個別電気機器の稼働状態x_t^-とt期の個別電気機器の規格電力P_t ^kの積和による、t期の総使用電力の予測値H_tx_t^^-と観測値との誤差y_t-H_tx_t^^-をベースに、カルマンゲインによって、状態変数x_t^^-を改定することで、t期までの情報をもちいた状態変数x_tの推定値x_t^⁺を次で求める。
これを、行列形式で表現すれば、以下となる。
最後に、状態変数の推定量の最小共分散行列Pの推定値を、t-1期までの情報をもちいた推定値P_t ^-からt期までの情報をもちいた推定値P_t ⁺に、次の式で更新する。

以上のプロセスにより、期のアウトプットとして、x_t^⁺,P_t ⁺およびH_tを出力し、t+1期に移る。これが、カルマンフィルタの第期の計算プロセスである。
図７は、これまで説明してきた、本実施例実施例の、NIALMのためのリアルタイム逆推定アルゴリズムでの全体の計算プロセスを、フローチャートにまとめた図である。

ここで、リアルタイム逆推定アルゴリズムの流れを再度まとめたものを、図７に従って以下に概説する。
１） 初期値を与える。初期値には、個別電気機器の稼働確率x₀^⁺、および、その共分散P₀ ⁺、また、状態方程式の誤差共分散Q_t, t=0,1,…、および、観測方程式の誤差共分散R_t, t=0,1,2,…も与える。状態方程式および観測方程式の誤差共分散は、すべての期で外生的に与える。

２） k-ミーンズアルゴリズムによって、総使用電力のジャンプがこれまでに存在している電気機器の規格電力かどうかを判断し、新たなジャンプであれば、新規の電気機器として機種リストに加え、その規格電力を新しい機種の規格電力として仮登録する。

３） 仮登録された個別電気機器の機種と規格電力を含め、機種リストに登録された電気機器の規格電力およびその分散を再計算する。既存の機種間、あるいは、新規の機種と既存の機種が、規格電力の再計算の結果、同一のものとみなすことができれば、当該の電気機器を機種リストから取り除く。これによって、個別電気機器の機種リスト、その規格電力リストを最新のリストに改定する。

４） k-ミーンズアルゴリズムの計算結果をもとに、総使用電力のジャンプが、最新の機種リスト、規格電力リストに登録された、どの規格電力をもつ、どの個別電気機器の稼働状態の変化に起因するジャンプなのか、の起こりやすさの確率を計算する。

５） 総使用電力のジャンプが、プラスか、マイナスか、0かによって、状態推移が異なるマルコフスイッチングモデルを構成し、総使用電力のジャンプの符号と、t-1期の個別電気機器の稼働状態の違いによって、t期への状態推移が異なる状態方程式を構成する。

６） t-1期までの情報をもちいた個別電気機器の稼働確率と、t期での個別電気機器の規格電力の推定値の情報を用いてt期の総使用電力を予測する。その予測値とt期の総使用電力の観測値との誤差にもとづいて、カルマンフィルタによって、状態変数と状態変数の共分散行列の最適推定の再帰的な更新を行い、t期までの情報をもちいた、t期の個別電気機器の稼働確率xt^+、および、その共分散行列Pt+を出力する。

本実施例のリアルタイム逆推定アルゴリズムは、次の特徴と意義を持っている。
1) 層使用電力の時間変動の１階階差、すなわち、ジャンプにのみ着目する点。
2) 電気機器の機種を企画電力の違いで識別する単純化を行っている点。
3) ジャンプを、クラスタの数が変化する可変的k−ミーンズモデルへの入力と捉え、EMアルゴリズムの解法を用いることで、１分間隔のリアルタイムでの機種の同定と推定を可能とした点。
4) 可変的なk−ミーンズモデルと状態空間モデリングのカルマンフィルタとを接合し、同定した個別電気機器の稼働状態の変化を、リアルタイムの動的推定可能としている点。
5) k−ミーンズモデルとカルマンフィルタの統合には、マルコフスイッチングの考え方を梃子としている点。
このように、本実施例の逆推定アルゴリズムは、リアルタイムの動的推定を可能にする実装を意図した実用性の高い、これまでにないアルゴリズムである。

次に、本実施例のリアルタイム逆推定アルゴリズムを、実際の過程環境データへ適用し、その評価を行った。
実際の過程環境データへの適用は、以下のように行った。
実際に生活している環境（東京都内）にNIALMを８日間（2013年3月12日〜19日）設置して、総使用電力を観測し、個別電気機器の種類と、規格電力を推定できるかを行った。NIALMは疑似家庭環境と同機を使用しており、１秒おきにデータを収集した。実験に協力していただいた家庭は、４人家族で、うち２人は子供である。実験期間の天候は図８に示す表の通りである。
本実施例では、今後、一般家庭に普及することを想定して、８日間を１分間隔のデータ11520個のデータを使用した。

この評価結果を図９〜１６に示す。
図９は、総使用電力の1分間のジャンプの規模別度数分布を表している。同図から明らかなように、観測されたジャンプは0を中心としてプラスとマイナスに対称的に描かれていることがわかる。それゆえ、総使用電力のジャンプはその機器がオン、オフになった可能性が非常に高く、そのジャンプが規格電力である可能性が高いといえる。
本実施例で構築した逆推定アルゴリズムの適用対象としては、実家庭環境での総使用電力の1分間隔の時間変動データを想定している。とくに、総使用電力の1分時間間隔の時間変動データの1階階差、ジャンプに着目し、それを逆推定アルゴリズムの入力ストリームとしている。
それは、個別電気機器がスイッチオンしたときに、規格電力のポジティブのジャンプが起こり、いずれかの時点で、個別電気機器がスイッチオフとなったときに、ネガティブのジャンプが起こり、個別電気機器の使用電力は0となる、との考え方である。これは、特許文献１に記載の技術のように総使用電力の波形にまで踏み込まずに、オンとオフのジャンプの動きにのみ着目する方法である。

図１０は、総使用電力規模別総電力使用量に占めるシェアを表している。この値の合計を求めると、総電力使用量8,941,000Wmが得られる。これを1ヶ月あたりの総電力使用量に換算すると、558.8kWh/月＝8941kWm÷60分÷8日×30日となる。これは、1世帯あたりの電力使用量544kWh/月（参照：（2013年3月）総務省家計調査）の値に近く、この被験者の家庭は平均的な電気の使い方をしており、NIALMによる電力推定の実験対象に適していることがわかる。

図１１は、８日間実環境データのリアルタイム処理アルゴリズムの結果で、使用開始から1200分までにかけてのものである。同図において横軸には時間を、また縦軸には総使用電力を示す。
なお、このアルゴリズムでは、状態変数の引き渡しについては行っていない。
図１２は、図１１のうち使用開始から150分までを切り出して拡大したものである。また、図１３は、図１１のうち150分から300分までを切り出して拡大したものである。
図１２から、使用開始から150分にかけては、総使用電力の観測値と推定値で誤差が生じている箇所がいくつかあるが、その後の150分以降は、図１３のように、観測値と推定値の間に大きな誤差は生じておらず、推定精度がよいことが分かった。

図１４は、個別電気機器を手動でオンとオフした疑似家庭環境での４時間の１分間隔のデータを用いて、本実施例を適用し、個別電気機器の種類および規格電力を推定し、総使用電力と個別電気機器の稼働状態を推定したものである。本実施例のカルマンフィルタでは、状態変数ベクトルの次元を可変としているが、状態変数が変化する場合、変化以前の状態変数を変化後にどのように引き継ぐかを考慮する必要がある。

図１４は、個別電気機器の種類の数が変わったときの状態変数の引き渡しを考慮しておらず、個別電気機器の稼働確率をゼロクリアしている。使用開始から150分にかけては、個別電気機器の種類や規格電力が未知であり、個別電気機器の種類が変化するにもかかわらず、状態変数の引き渡しが考慮されていないため、総使用電力の観測値と推定値とが大きく異なる時点がいくつもある。
ただし、150分以降は個別電気機器の種類の数が変わったときの状態変数の引き渡しを考慮しなくても、個別電気機器の種類の数は学習で安定した状態に達し、総使用電力の観測値と推定値との間に大きな違いは生じていない。

これに対し、図１５は、個別電気機器の種類および規格電力を推定し、個別電気機器の種類の数が変わったときの状態変数の引き渡しを考慮した上で、総使用電力と個別電気機器の稼働状態を推定したものである。
図１３の状態変数の引き渡しを考慮していない場合と比べると、使用開始から、総使用電力の観測値と推定値との間に大きな違いが見られず、初期の段階から予測精度の向上が図られている。

以上、説明したように、実施例１の個別電気機器稼働状態推定装置およびその方法は、以下の効果を有する。
すなわち、実施例１の個別電気機器稼働状態推定装置およびその方法にあっては、個々の電気機器にセンサを設けることなく、また事前にどのような種類の電気機器が何台あるかといった情報がなくても、個々の電気機器の使用状況（どの電気機器の機種があり、どの電気機器がどのような稼働状態にあるか）をリアルタイムでモニタリングを行うことができる。したがって、安価になり、一般家庭等への適用も容易となって、その結果、電気機器の使用状況の改善や電力供給制御の改善等を実施して省エネルギ化を図ることが可能となる。

また、電気機器機種特定部6aで行うクラスタリングにはk-ミーンズアルゴリズムを用い、このk-ミーンズアルゴリズムで用いる前記個別電気機器の機種に対応するクラスタの数kを可変としたので、途中で新しい電気機器を追加したり、あるいは途中でそれまで用いていた電気機器を廃棄したりして、モニタリングの対象となる電気機器の機種が異なったとしても、その都度、機種リストのデータの追加入力や削除を行う必要がなくなり、その工数も不要となる。この結果、素人でも、本実施例の個別電気機器稼働状態推定装置を利用することが可能となる。

また、k-ミーンズアルゴリズムが、期待値ステップおよび最大化ステップを有するEMアルゴリズムを用いて、入力された前記ジャンプ電力に基づいて未知のk個の正規分布からの実現値であるときのk個の正規分布の平均と分散を推定し、前記入力されたジャンプ電力がどの正規分布からの実現値であるかを判定するようにしたので、電気機器の増減にかかわらず、事前情報を入力しなくても、個別電気機器について必要な情報を入手することができる。

また、EMアルゴリズムでの期待値ステップおよび最大化ステップの、時間変動の時間間隔ごとの繰り返し計算は、通常何回も繰り返し計算を行うが、本実施例では１ステップの再計算のみとしたので、EMアルゴリズムの尤度の単調性から、時間経過とともに推定精度が上がり、推定精度を十分維持しながらも、計算コストを低く抑えることが可能となる。

また、通常、k-ミーンズアルゴリズムでは、すべてのデータを用いるが、本実施例の最大化ステップでは、用いる過去のデータの長さ、逆戻る長さは、これらより前から得られた過去のデータより少ない長さに指定するようにしたので、それまで得た全データを用いることなく、計算が可能となり、計算コストを低く抑えることが可能となる。

また、EMアルゴリズムでは、期待値ステップでの期待値の算出に、この期待値の事後確率が期待値の事前分布と尤度との積に比例するとしたベイジアン学習に基づく事前分布を導入し、総使用電力とジャンプ電力の過去のデータから得られた、個別電気機器のスイッチオン、オフの発生確率の時間変動、使用環境の違いによる変動、使用日変動のうちの少なくとも一つの変動のパターンの学習結果の推定に利用するようにした、できるようにしたので、最尤推定を行うことが可能となる。

また、各電気機器の稼働確率の状態変数に対応させるマルコフスイッチングモデルと、状態方程式の標準形に、前記総使用電力に対して前記ジャンプ電力を対応させた各入力変数を持たせ、前記各電気機器の稼働確率とした状態変数を用いるカルマンフィルタと、を備えているので、k-ミーンズでの推定結果を、マルコフスイッチングモデルを用いてカルマンフィルタで利用できるように変換することができる。その結果、個別電気機器の稼働状態を動的にリアルタイムで推測することが可能となる。

また、カルマンフィルタの状態変数である状態変数ベクトルの次元は、通常不変であるが、本実施例でのカルマンフィルタの状態変数である状態変数ベクトルの次元は、可変としたので、途中で新しい電気機器を追加したり、あるいは途中でそれまで用いていた電気機器を廃棄したりし場合でも、その都度、機種リストのデータの追加入力や削除を行う必要がなくなり、その工数も不要となる。

以上、本発明を上記実施例に基づき説明してきたが、本発明は上記実施例に限られず、本発明の要旨を逸脱しない範囲で設計変更等があった場合でも、本発明に含まれる。

たとえば、本発明の個別電気機器稼働状態推定装置およびその方法は、上記実施例では１軒の一般家庭について適用したが、一般家庭でなくてもよく、あるいは集団住宅をまとめたものに適用するようにしてもよい。
あるいは、本発明の個別電気機器稼働状態推定装置およびその方法によるモニタリングの対象としては、家屋１軒を複数のエリアに分けて、それぞれのエリアに分電盤などを設け、個別あるいは統合して適用するようにしてもよい。

また、電力計は通常の家庭に取り付けられているものに限られず、その他のものを利用するようにしてもよい。また、各プロセスでの演算は、上記演算式に限られない。

また、上記実施例では、ジャンプ電力を総使用電力の1階階差から求めたが、発明ではこれに限られない。

本発明の個別電気機器稼働状態推定装置およびその方法は、モニタリングの対象となる複数の電気機器について、その稼働状態を推測する種々のシステムに利用することが可能となる。

1 家屋
2 配線
3a〜3n 電気機器
4 供給電線
5 電力計（総使用電力測定手段）
6 稼働状態推定部
6a 時間変動算出部（時間変動算出手段）
6b 稼働電気機器機種特定部（稼働電気機器機種特定部）
6c 対応尤度推定部（対応尤度推定手段）
6d 個別電気機器稼働状態推定部（個別電気機器稼働状態推定手段）
6e マルコフスイッチングモデル
6f カルマンフィルタ
7 ディスプレイ
8 分電盤（総使用電力測定手段）

本発明は、複数の電気機器のうち、稼働している電気機器の種類（機種）と、その稼働状態とを同時にリアルタイムでモニタリング可能な個別電気機器稼働状態推定装置、およびその方法に関する。

すなわち、本発明の個別電気機器稼働状態推定装置は、
複数の電気機器の総使用電力を測定する総使用電力測定手段と、
総使用電力測定手段で測定した総使用電力に基づき、この総使用電力の時間的変動であるジャンプ電力を算出する時間変動算出手段と、
時間変動算出手段で算出したジャンプ電力に基づいて新たな電気機器の機種として既存の機種リストに加えるか否かを判定し、新たな電気機器の機種として既存の機種リストに加える場合には個別電気機器の機種別の規格電力の再計算を行い、必要に応じて既存の機種リストから余分な機種を削除し、新たな電気機器の機種を加えない場合には既存の機種リストを維持する、前記電気機器の機種のクラスタリングを行う電気機器機種特定手段と、
時間変動算出手段で算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定する対応尤度推定手段と、
対応尤度推定手段で得られた、現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの個別電気機器の稼働状態の変化を推定して、現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と総使用電力の測定データから個別電気機器の機種の推定、規格電力の推定、および稼働状態の推定を最適化し、測定された総使用電力の時間変動の時間間隔で刻々変化する個別電気機器の稼働確率を動的に推定する個別電気機器稼働状態推定手段と、
を備えたことを特徴とする。

一方、本発明の個別電気機器稼働状態推定方法は、
複数の電気機器の総使用電力を測定するステップと、
測定した総使用電力に基づき、この総使用電力の時間的変動であるジャンプ電力を算出するステップと、
算出したジャンプ電力に基づいて新たな電気機器の機種として既存の機種リストに加えるか否かを判定し、新たな電気機器の機種として既存の機種リストに加える場合には個別電気機器の機種別の規格電力の再計算を行い、必要に応じて既存の機種リストから余分な機種を削除し、新たな電気機器の機種を加えない場合には既存の機種リストを維持する、前記電気機器の機種のクラスタリングを行うステップと、
算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器の機種のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定するステップと、
現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの個別電気機器の稼働状態の変化を推定して、現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と総使用電力の測定データから個別電気機器の機種の推定、規格電力の推定、および稼働状態の推定を最適化し、測定された総使用電力の時間変動の時間間隔で刻々変化する個別電気機器の稼働確率を動的に推定するステップと、
を有することを特徴とする。

本発明の個別電気機器稼働状態推定装置およびその方法にあっては、個々の電気機器にセンサを設けることなく、また事前にどのような種類の電気機器が何台あるかといった情報がなくても、モニタリングの対象となる電気機器の機種ごとの使用状況をリアルタイムでモニタすることができる。

ステップS10では、個別電気機器稼働状態推定部6dが、マルコフスイッチングモデル6eを用いて、ステップS8で得られた、現在の個別電気機器の稼働状態のもとで総使用電力の変動ジャンプから得られる個別電気機器のスイッチオン/スイッチオフの事象生起に関する尤度が入力されたとき、どのように個別電気機器の稼働状態が変化するか、を推定し、現在の個別電気機器の稼動状態の推定を更新する。この推定結果は、カルマンフィルタ6fへ出力する。続いて、ステップS11へ進む。

Claims (16)

1. 複数の電気機器の総使用電力を測定する総使用電力測定手段と、
   該総使用電力測定手段で測定し総使用電力に基づき、該総使用電力の時間的変動であるジャンプ電力を算出する時間変動算出手段と、
   該時間変動算出手段で算出したジャンプ電力に基づいて新たな電気機器として既存の機種リストに加えるか否かを判定し、新たな電気機器として前記既存の機種リストに加える場合には個別電気機器別の規格電力の再計算を行い、必要に応じて前記既存の機種リストから余分な機種を削除し、新たな電気機器を加えない場合には前記既存の機種リストを維持するクラスタリングを行う電気機器機種特定手段と、
   前記時間変動算出手段で算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定する対応尤度推定手段と、
   該対応尤度推定手段で得られた、現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、前記個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの前記個別電気機器の稼働状態の変化を推定して、前記現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と前記総使用電力の測定データから前記個別電気機器の機種の推定、前記規格電力の推定、および前記稼働状態の推定を最適化し、測定された前記総使用電力の時間変動の時間間隔で刻々変化する前記個別電気機器の稼働確率を動的に推定する個別電気機器稼働状態推定手段と、
   を備えたことを特徴とする個別電気機器稼働状態推定装置。
2. 請求項１に記載の個別電気機器稼働状態推定装置において、
   前記電気機器機種特定手段で行うクラスタリングにはk-ミーンズアルゴリズムを用い、該k-ミーンズアルゴリズムで用いる前記個別電気機器の機種に対応するクラスタの数kを可変とした、
   ことを特徴とする個別電気機器稼働状態推定装置。
3. 請求項２に記載の個別電気機器稼働状態推定装置において、
   前記k-ミーンズアルゴリズムは、期待値ステップおよび最大化ステップを有するEMアルゴリズムを用いて、入力された前記ジャンプ電力に基づいて未知のk個の正規分布からの実現値であるときのk個の正規分布の平均と分散を推定し、前記入力されたジャンプ電力がどの正規分布からの実現値であるかを判定する、
   ことを特徴とする個別電気機器稼働状態推定装置。
4. 請求項３に記載の個別電気機器稼働状態推定装置において、
   前記EMアルゴリズムでの期待値ステップおよび最大化ステップの、時間変動間隔ごとの繰り返し計算は、１ステップの再計算のみとした、
   ことを特徴とする個別電気機器稼働状態推定装置。
5. 請求項３又は４に記載の個別電気機器稼働状態推定装置において、
   前記最大化ステップに用いる、平均および分散の過去のデータの長さ、逆戻る長さは、該過去のデータより前から得てきた過去のデータより少ない長さに指定する、
   ことを特徴とする個別電気機器稼働状態推定装置。
6. 請求項３乃至５のうちの１項に記載された個別電気機器稼働状態推定装置において、
   前記EMアルゴリズムは、前記期待値ステップでの期待値の算出に、該期待値の事後確率が前記期待値の事前分布と尤度との積に比例するとしたベイジアン学習に基づく前記事前分布を導入し、前記総使用電力と前記ジャンプ電力の過去のデータから得られた、前記個別電気機器のスイッチオン、オフの発生確率の時間変動、使用環境の違いによる変動、使用日変動のうちの少なくとも一つの変動のパターンの学習結果の推定に利用するようにした、
   ことを特徴とする個別電気機器稼働状態推定装置。
7. 請求項２乃至６のうちのいずれか１項に記載の個別電気機器稼働状態推定装置において、
   各電気機器の稼働確率の状態変数に対応させるマルコフスイッチングモデルと、
   状態方程式の標準形に、前記総使用電力に対して前記ジャンプ電量を対応させた各入力変数を持たせ、前記各電気機器の稼働確率とした状態変数を用いるカルマンフィルタと、を備えている、
   ことを特徴とする個別電気機器稼働状態推定装置。
8. 請求項７に記載の個別電気機器稼働状態推定装置において、
   前記カルマンフィルタの状態変数を、次元が可変である状態変数ベクトルとした、
   ことを特徴とする個別電気機器稼働状態推定装置。
9. 複数の電気機器の総使用電力を測定するステップと、
   前記測定した総使用電力に基づき、該総使用電力の時間的変動であるジャンプ電力を算出するステップと、
   前記算出したジャンプ電力に基づいて新たな電気機器として既存の機種リストに加えるか否かを判定し、新たな電気機器として前記既存の機種リストに加える場合には個別電気機器別の規格電力の再計算を行い、必要に応じて前記既存の機種リストから余分な機種を削除し、新たな電気機器を加えない場合には前記既存の機種リストを維持するクラスタリングを行うステップと、
   前記算出したジャンプ電力に基づいてアップおよびダウンのジャンプがどの電気機器のスイッチオン及びスイッチオフの事象の生起に対応しているかの尤度を推定するステップと、
   現在の個別電気機器の稼働状態のもとで総使用電力のジャンプから得られる、前記個別電気機器のスイッチオン及びスイッチオフの事象生起に関する尤度が入力されたときの該個別電気機器の稼働状態の変化を推定して、前記現在の個別電気機器の稼働状態の推定を更新し、この更新された個別電気機器の稼働状態の推定結果と規格電力の推定結果を用いて予測される総使用電力の予測値と前記総使用電力の測定データから前記個別電気機器の機種の推定、前記規格電力の推定、および前記稼働状態の推定を最適化し、測定された前記総使用電力の時間変動の時間間隔で刻々変化する前記個別電気機器の稼働確率を動的に推定するステップと、
   を有する別電気機器稼働状態推定方法。
10. 請求項９に記載の個別電気機器稼働状態推定方法において、
    前記クラスタリングは、k-ミーンズアルゴリズムを用い、該k-ミーンズアルゴリズムで用いる前記個別電気機器の機種に対応するクラスタの数kを可変とした、
    ことを特徴とする個別電気機器稼働状態推定方法。
11. 請求項１０に記載の個別電気機器稼働状態推定方法において、
    前記k-ミーンズアルゴリズムは、期待値ステップおよび最大化ステップを有するEMアルゴリズムを用いて、入力された前記ジャンプ電力に基づいて未知のk個の正規分布からの実現値であるときのk個の正規分布の平均と分散を推定し、前記入力されたジャンプ電力がどの正規分布からの実現値であるかを判定する,
    ことを特徴とする個別電気機器稼働状態推定方法。
12. 請求項１１に記載の個別電気機器稼働状態推定方法において、
    前記EMアルゴリズムでの期待値ステップおよび最大化ステップの、時間変動間隔ごとの繰り返し計算は、１ステップの再計算のみとした、
    ことを特徴とする個別電気機器稼働状態推定方法。
13. 請求項１１又は１２に記載の個別電気機器稼働状態推定方法において、
    前記最大化ステップに用いる、平均および分散の過去のデータの長さ、逆戻る長さは、該過去のデータより前から得てきた過去のデータ過去のデータより少ない長さに指定する、
    ことを特徴とする個別電気機器稼働状態推定方法。
14. 請求項１１乃至１３のうちの１項に記載された個別電気機器稼働状態推定方法において、
    前記EMアルゴリズムは、前記期待値ステップでの期待値の算出に、該期待値の事後確率が前記期待値の事前分布と尤度との積に比例するとしたベイジアン学習に基づく前記事前分布を導入し、前記総使用電力と前記ジャンプ電力の過去のデータから得られた、前記個別電気機器のスイッチオン、オフの発生確率の時間変動、使用環境の違いによる変動、使用日変動のうちの少なくとも一つの変動のパターンの学習結果の推定に利用するようにした、
    ことを特徴とする個別電気機器稼働状態推定方法。
15. 請求項９乃至１４のうちのいずれか１項に記載の個別電気機器稼働状態推定方法において、
    マルコフスイッチングモデルを用いて各電気機器の稼働確率の状態変数に対応させるステップと、
    カルマンフィルタの状態方程式の標準形に、前記総使用電力に対して前記ジャンプ電量を対応させた各入力変数を持たせ、前記カルマンフィルタで前記各電気機器の稼働確率とした状態変数を用いると、
    を有することを特徴とする個別電気機器稼働状態推定方法。
16. 請求項１５に記載の個別電気機器稼働状態推定方法において、
    前記カルマンフィルタの状態変数を、次元が可変である状態変数ベクトルとした、
    ことを特徴とする個別電気機器稼働状態推定方法。