(57)【特許請求の範囲】
【請求項１】
複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換するデータ取得手段と、
非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求めるパラメータ推定手段と
を備え、
前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求める
データ処理装置。
【請求項２】
前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルの各状態が取り得る状態数を２とし、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータを求める
請求項１に記載のデータ処理装置。
【請求項３】
前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルのファクタmの波形のパターンに対応する観測確率のパラメータW(m)の、前記確率モデルの各状態が取り得る状態数に対応する複数列のいずれか一列が常に０であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求める
請求項１または２のいずれかに記載のデータ処理装置。
【請求項４】
前記確率モデルの前記ファクタ１乃至ファクタMのうちの所定の複数のファクタが同一の前記波形に対応するものであることを判定する同一判定手段をさらに備え、
前記同一判定手段は、前記確率モデルの各ファクタの全組み合わせのうち、必要条件を常に満たす組み合わせを同一候補として選別し、選別された同一候補のうち、ヒント条件を所定回数以上満たしたことがある組み合わせを同一と判定することにより、所定の複数のファクタが同一の前記波形に対応するものであることを判定する
請求項１乃至３のいずれかに記載のデータ処理装置。
【請求項５】
現在時刻から所定時間経過後の未来の前記ファクタの状態を予測する状態予測手段をさらに備える
請求項１乃至４のいずれかに記載のデータ処理装置。
【請求項６】
前記ファクタは第１の機器に対応するものであり、前記状態予測手段による未来の前記ファクタの状態の予測結果に基づいて、前記第１の機器に関連する第２の機器を制御する機器制御手段をさらに備える
請求項５に記載のデータ処理装置。
【請求項７】
取得された前記時系列データに基づいて推定された前記ファクタの観測確率のパラメータW(m)を、予め記憶されている所定の製品種類の電流波形パターンと比較することにより、前記複数の波形それぞれに対応する複数の機器を特定する機器特定手段をさらに備える
請求項１乃至６のいずれかに記載のデータ処理装置。
【請求項８】
前記ファクタに対応する機器をユーザに指定させることにより、前記複数の波形それぞれに対応する複数の機器を特定する機器特定手段をさらに備える
請求項１乃至６のいずれかに記載のデータ処理装置。
【請求項９】
前記パラメータ推定手段の結果に基づいて、前記機器特定手段により特定された前記複数の機器の状態を表示する表示手段をさらに備える
請求項７または８に記載のデータ処理装置。
【請求項１０】
前記データ取得手段は、他の測定値の位相に同期して、所定の周期で前記複数の波形の合計値の時系列データを取得する
請求項１乃至９のいずれかに記載のデータ処理装置。
【請求項１１】
データ処理装置が、
複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換し、
非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求め、
EMアルゴリズムによるパラメータ推定処理では、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求める
ステップを含むデータ処理方法。
【請求項１２】
コンピュータを、
複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換するデータ取得手段と、
非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求めるパラメータ推定手段として機能させ、
前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求める
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、データ処理装置、データ処理方法、およびプログラムに関し、特に、取得される電流情報から電気機器を推定する場合において、パラメータの算出及び調整が簡単であり、かつ、事前のデータベースが不要な手法を確立することができるようにするデータ処理装置、データ処理方法、およびプログラムに関する。
【背景技術】
【０００２】
分電盤で計測した電流の情報から、その先につながっている電気機器を推定する技術は、 Non-intrusive load monitoring (以下、NILMという。)と呼ばれ、1980年代から研究されている。NILMは、個々の電気機器（負荷）それぞれに測定機器を必要とせず、一箇所の計測結果だけに基づいて、そこから先につながっている電気機器全ての状態を把握できることが大きな利点となっている。
【０００３】
NILMの代表的なものとして、例えば、特許文献１には、電流と電圧の計測から、有効電力（real power）及び無効電力（reactive power）を算出し、その変化量をクラスタリングすることで、電気機器の同定を行う技術が開示されている。変化量を取るのは、電気機器のオン、オフ時に、計測される有効電力及び無効電力が変化するからである。
【０００４】
図１は、特許文献１において図８として示されている図である。図１では、有効電力及び無効電力を軸とする２次元平面に、冷蔵庫（Refrigerator）と暖房（Heater）のオン、オフ時の有効電力と無効電力がプロットされている。図１によれば、電気機器のオンとオフが点対称の位置にプロットされていることがわかる。
【０００５】
特許文献１の方法は、オン、オフ時の差分を取るので、変化した瞬間の情報しか利用していない。また、変化点の検出器（Change detector）が必要となるが、この変化点検出器が失敗する（オン若しくはオフを取り逃す、または、過剰に変化を誤検出する）と、後段の処理全体が失敗してしまう。
【０００６】
即ち、特許文献１の方法には、次のような問題がある。第１に、オン、オフ時の差分を取るので、変化した瞬間の情報しか利用していない。第２に、変化点検出のしきい値調整が難しく、変化点検出器（net change detection）が失敗すると後段の処理が全て失敗してしまう。第３に、８０年代の当時の家電はシンプルな負荷であるものが多かったため適用できたが、現代の電気機器ではオン、オフだけに分類できず、うまく機能しないものが多い。
【０００７】
近年の複雑な電力消費を行う電気機器に対応するため、NILMの側も何らかの複雑な処理を行う必要が出てきている。その試みとして、識別モデル（判別モデル、Classification）を用いる方法が多く提案されている。識別モデルにSupport Vector MachineなどのLMC（Large Margin Classifier）を利用した方法として、例えば、特許文献２及び３に開示された方法がある。
【０００８】
AdaBoostやSupport Vector Machineなどの識別モデルは、特徴量をうまく選択し、かつ学習用のサンプルデータが十分に多い場合は、非常に高い識別性能を示すことが知られているため、この手法は精度向上には有効であるとも思われる。しかし、一方で、識別モデルはHMMなどの生成モデルとは違い、事前に学習データを用意して学習を済ませておく必要があり、更にその学習結果をデータベースとして保持しておく必要がある。つまり、未知の電気機器に対応することはできないという欠点がある。
【０００９】
特許文献４及び非特許文献１のように、単純な線形モデルで判別を行うものもあるが、事前のデータベースが必要であるという点では、特許文献２及び３に開示された技術と同じ問題を抱えている。事前に用意したデータベースを用いる手法は、その他、特許文献５や６として提案されているものがある。
【００１０】
結局のところ、上述した従来の手法は、精度を取るか、または、事前登録が不要な手法を取るか、のトレードオフになっている。近年、家庭の電気機器は非常に多様化しており、事前学習が必要である識別モデルでは、事実上、家庭での使用に耐えられないと考えられるため、やはり、事前登録不要な手法の方が望ましい。
【００１１】
そこで、事前学習が必要な識別モデルではなく、生成モデルを使う試みも既に行われている。例えば、生成モデルとしてHidden Markovモデル（HMM）を適用したものが提案されている（例えば、非特許文献２，３参照）。
【００１２】
しかし、生成モデルとして単純なHMMを適用した場合には、電気機器の数が増えると状態数が爆発し（膨大となり）、現実的なシステムを構築することができないという問題がある。例えば、個々の電気機器が、仮にオンとオフの２状態をもち、電気機器の数をnとしたとき、必要な状態数は、２ｎである。更に、状態遷移確率のサイズはその二乗（２ｎ）２となる。一般家庭の全ての電気機器が仮に２０個あったとして（近年では決して多いとは言えない）、必要な状態数は２20=1,048,576 であり、状態遷移確率は1,099,511,627,776 の大きさとなる。これはテラオーダーのサイズであり、近年のパーソナルコンピュータでも取り扱うのは困難である。
【００１３】
ちなみに、特許文献１の方法も、クラスタリングをベースとしており、原始的な生成モデルと考えることができるため、事前登録は不要である。特許文献１の方法のような、確率モデルによるアプローチが一般化する以前の時代の、モデル化をせずに問題を解く方法は、ヒューリスティックな方法と呼ばれる。ヒューリスティックな方法は、ファーストステップとしては良いが、方法を拡張していくうちに閾値などのパラメータがどんどん増えていき、それらの調整が困難になるという問題を抱えている。
【００１４】
近年、コンピュータを利用した自動認識技術は、確率モデルの導入により様々な難しい問題が解決できるようになってきている。確率モデルによりうまくモデル化できれば、ほとんどのパラメータが最尤推定（ML推定、Maximum Likelihood）や事後確率最大化（MAP推定、Maximum A Posteriori）、最小識別誤差（MCE、Minimum Classification Error）などにより求めることが可能である。Support Vector Machineのような識別モデルや、HMMのような生成モデルを利用することが、確率モデルによるモデル化に相当する。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】米国特許第４８５８１４１号明細書
【特許文献２】特開２００１−３３０６３０号公報
【特許文献３】国際公開第０１／０７７６９６号パンフレット
【特許文献４】特開２００８−０３９４９２号公報
【特許文献５】特開２００６−０１７４５６号公報
【特許文献６】特開２００９−２５７９５２号公報
【非特許文献】
【００１６】
【非特許文献１】稲垣伸吉、江上司、鈴木達也、中村久栄、伊藤公一、”電気機器の非侵入型稼働状態モニタリングシステム  ―稼働の離散状態に着目した整数計画法による解法―”、計測自動制御学会 離散事象システム研究会講演論文集 第42回、pp.33-38、2007年
【非特許文献２】Bons M., Deville Y., Schang D. 1994. Non-intrusive electrical load monitoring using Hidden Markov Models.Third international Energy Efficiency and DSM Conference, October 31, Vancouver, Canada. , p. 7
【非特許文献３】中村久栄、伊藤公一、鈴木達也、「隠れマルコフモデルに基づいた電気機器の稼働状況モニタリングシステム」、電気学会論文誌Ｂ、Vol. 126、No. 12、pp. 12231229、2006年
【発明の概要】
【発明が解決しようとする課題】
【００１７】
以上のことから、次の３点を解決するNILMの手法が望まれる。第１に、ヒューリスティックな方法が抱えるようなパラメータ調整が困難であるという問題を解決すること、即ち、パラメータ調整が簡単であることが望まれる。第２に、識別モデルのように事前のデータベースが必要であると、近年益々増大する新しい電気機器に対応することが困難になることから、事前のデータベースが不要であることが望まれる。第３に、状態数が爆発してしまうパラメータ推定アルゴリズムでは実質的に解くことができないので、状態数がパラメータ算出可能な現実的な個数であることが必要である。
【００１８】
本発明は、このような状況に鑑みてなされたものであり、取得される電流情報から電気機器を推定する場合において、パラメータの算出及び調整が簡単であり、かつ、事前のデータベースが不要な手法を確立することができるようにするものである。
【課題を解決するための手段】
【００１９】
本発明の一側面のデータ処理装置は、複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換するデータ取得手段と、非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求めるパラメータ推定手段とを備え、前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求める。
【００２０】
本発明の一側面のデータ処理方法は、データ処理装置が、複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換し、非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求め、EMアルゴリズムによるパラメータ推定処理では、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求めるステップを含む。
【００２１】
本発明の一側面のプログラムは、コンピュータを、複数の波形の合計値の時系列データを取得し、取得した前記時系列データを、非負のデータに変換するデータ取得手段と、非負のデータに変換された前記時系列データに基づいて、前記複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータを求めるパラメータ推定手段として機能させ、前記パラメータ推定手段は、EMアルゴリズムによるパラメータ推定処理において、前記確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、前記確率モデルが、前記時系列データが表す前記波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、前記モデルパラメータとしての観測確率のパラメータW(m)を求めるものである。
【００２２】
本発明の一側面においては、複数の波形の合計値の時系列データが取得され、取得された時系列データが、非負のデータに変換され、非負のデータに変換された時系列データに基づいて、複数の波形それぞれが示す状態を確率モデルによりモデル化したときのモデルパラメータが求められる。EMアルゴリズムによるパラメータ推定処理では、確率モデルのファクタm（m＝１乃至Mの正の整数）の波形のパターンに対応する観測確率のパラメータW(m)が非負であるという制約条件の下で、確率モデルが、時系列データが表す波形の合計値のパターンを説明する度合いである尤度関数を最大化することにより、モデルパラメータとしての観測確率のパラメータW(m)が求められる。
【発明の効果】
【００２３】
本発明の一側面によれば、取得される電流情報から電気機器を推定する場合において、パラメータの算出及び調整が簡単であり、かつ、事前のデータベースが不要な手法を確立することができる。
【図面の簡単な説明】
【００２４】
【図１】従来の特許文献１の方法を説明する図である。
【図２】本発明を適用したデータ処理装置としての電気機器推定装置の概要について説明する。
【図３】通常のHMMとFactorial HMMとの違いを説明する図である。
【図４】Factorial HMMと図２の各電気機器との対応を示す図である。
【図５】図２の電気機器推定装置の構成例を示すブロック図である。
【図６】センサ部の詳細構成例を示すブロック図である。
【図７】電流測定部と電圧測定部の例を示す図である。
【図８】電流測定部により測定された電流値の時系列データの例を示す図である。
【図９】Factorial HMMのモデルパラメータとNILMにおける事象との対応関係を示す図である。
【図１０】完全に分解した変分法によるパラメータ推定処理を説明するフローチャートである。
【図１１】図１０のEステップ処理の詳細を説明するフローチャートである。
【図１２】図１０のMステップ処理の詳細を説明するフローチャートである。
【図１３】構造を残した変分法によるパラメータ推定処理を説明するフローチャートである。
【図１４】図１３のEステップ処理の詳細を説明するフローチャートである。
【図１５】パラメータ推定部によるMステップ処理について説明するフローチャートである。
【図１６】本来の状態数Kが３である電気機器のパラメータ推定結果を示す図である。
【図１７】同一機器判定部による同一機器判定処理を説明するフローチャートである。
【図１８】図５の電気機器推定装置全体の処理を説明するフローチャートである。
【図１９】本発明を適用したコンピュータの一実施の形態の構成例を示すブロック図である。
【発明を実施するための形態】
【００２５】
［本発明の概要］
初めに、図２を参照して、本発明を適用したデータ処理装置としての電気機器推定装置の概要について説明する。
【００２６】
住宅等においては、電力会社から供給される電気は、まず分電盤１１に引き込まれ、分電盤１１から、住宅内の各場所に設置された電気機器１２に供給される。例えば、図２の例では、分電盤１１から、電気機器１２としての、照明装置（電球）１２−１、エアコンディショナ１２−２、洗濯機１２−３、冷蔵庫１２−４、及びテレビジョン受像機１２−５に供給されている。
【００２７】
本発明を適用した電気機器推定装置１は、分電盤１１の大元の２次側において、住宅内の各場所に設置された複数の電気機器１２の使用状態の組み合わせでなる消費電流の合計値を測定することで、複数の電気機器１２の稼働状態を推定する。電気機器推定装置１は、推定結果に基づいて、各電気機器１２の現在の稼働状態を表示したり、現在より所定時間経過後の未来の電気機器１２の稼働状態を予測する。
【００２８】
［Factorial HMMについて］
電気機器推定装置１は、NILM の解析手段としてFactorial HMM（Hidden Markov Model）を用いて、各電気機器１２の稼働状態を推定する。換言すれば、電気機器推定装置１は、各電気機器１２の稼働状態を、Factorial HMMによりモデル化したモデルパラメータを求める。
【００２９】
そこで、初めに、Factorial HMMについて簡単に説明する。図３は、通常のHMMとFactorial HMMをグラフィカルモデルで表現した図である。
【００３０】
図３Aが、通常のHMMを表現したグラフィカルモデルであり、図３BがFactorial HMMを表現したグラフィカルモデルである。
【００３１】
通常のHMMでは、時刻tの観測データYtに対して、１つの状態変数Stが対応する。この通常のHMMとFactorial HMMが異なる点は、状態変数StがSt (１)，St (２)，St (３)，・・・St (m)，・・・St (M)と、複数（図３ではM個）存在し、それらの複数の状態変数St (１)乃至St (M)から、１つの観測データYtが生成されることである。
【００３２】
図４は、図３BのFactorial HMMを、図２に示した各電気機器１２と対応させて示した図である。
【００３３】
Factorial HMMのM個の状態変数S(１) 乃至S(M)のそれぞれが各電気機器１２に対応する。また、状態変数S (m)の状態値は、電気機器１２の状態（例えば、オン、オフの２状態）に対応する。
【００３４】
より具体的には、M個の状態変数S(１) 乃至S(M)のうち、１番目の状態変数S(１)の時間経過に応じた状態値S１ (１)乃至St (１)が、所定の電気機器１２（例えば、冷蔵庫１２−４）の状態に対応する。また、２番目の状態変数S(２)の時間経過に応じた状態値S１ (２)乃至St (２)が所定の電気機器１２（例えば、テレビジョン受像機１２−５）の状態に対応する。同様に、m番目の状態変数S(m)の時間経過に応じた状態値S１ (m)乃至St (m)が所定の電気機器１２（例えば、洗濯機１２−３）の状態に対応する。
【００３５】
また、住宅内の各場所に設置された複数の電気機器１２の使用状態の組み合わせでなる消費電流の合計値が、観測データY１乃至Ytとして得られる。
【００３６】
なお、以下では、M個の状態変数S(１) 乃至S(M)のうちの、m番目の状態変数S(m)を、m番目のファクタ、または、ファクタmとも記述する。
【００３７】
Factorial HMMの詳細は、Zoubin Ghahramani, and Michael I. Jordan, Factorial Hidden Markov Models’, Machine Learning Volume 29, Issue 2-3 ,Nov./Dec. 1997（以下、文献Xという。）に、記載されている。
【００３８】
［電気機器推定装置１の構成例］
図５は、電気機器推定装置１の構成例を示すブロック図である。
【００３９】
電気機器推定装置１は、センサ部２１、パラメータ推定部２２、データベース２３、同一機器判定部２４、状態予測部２５、機器特定部２６、及び表示部２７により構成される。
【００４０】
センサ部２１は、観測データYt（ｔ＝１，２，・・・，T）としての、住宅内の各場所に設置された複数の電気機器１２の使用状態の組み合わせでなる消費電流の合計値を測定（取得）し、パラメータ推定部２２に供給する。
【００４１】
パラメータ推定部２２は、複数の電気機器１２の消費電流の合計値の時系列データである観測データ｛Y１，Y２，Y３，・・・Yt，・・YT｝に基づいて、各電気機器１２の稼働状態をFactorial HMMによりモデル化したモデルパラメータを求める。Factorial HMMの学習処理により得られたモデルパラメータは、データベース２３に保存される。
【００４２】
また、パラメータ推定部２２は、センサ部２１から新たな観測データYtが供給された場合、モデルパラメータをデータベース２３から取得し、更新する。即ち、新たな観測データYtに基づいて、現在のモデルパラメータが更新される。
【００４３】
同一機器判定部２４は、M個のファクタのうち、同一の電気機器１２をモデル化した複数のファクタを検出し、検出結果をデータベース２３に記憶させる。換言すれば、同一機器判定部２４は、M個のファクタのうち、第１のファクタm１と第２のファクタm２（m１≠m２）が同一の電気機器１２を表現したものであるかを判定し、判定結果をデータベース２３に登録する。
【００４４】
Factorial HMM自体は、汎用の時系列データのモデル化手法であり、NILM以外の様々な問題に適用可能なものであるため、従来のFactorial HMMを用いた推定手法そのままでは、NILMにうまく適用できない。その問題の１つが、１つの電気機器１２が複数のファクタでモデル化される場合があることである。そこで、同一機器判定部２４は、１つの電気機器１２が複数のファクタで表現されている場合、その複数のファクタが同一の電気機器１２に対応するものであることを検出する。
【００４５】
状態予測部２５は、データベース２３に記憶されているモデルパラメータを用いて、現在時刻から所定時間経過後の未来のファクタm（の電気機器１２）の状態を予測する。Factorial HMMは、言うまでもなくHMMをベースとした確率モデルであるため、未来の時刻の状態確率を確率的に予測することができる。状態予測部２５は、現在時刻から所定時間経過後の未来のファクタmの状態の予測結果を、例えば、図示せぬ関連機器制御部等に出力する。これにより、関連機器制御部は、未来のファクタm（の電気機器１２）の状態の予測結果に基づいて、ファクタmの電気機器１２に関連する他の電気機器１２を制御することができる。
【００４６】
機器特定部２６は、ファクタmがどの電気機器１２に対応するのかを特定する。即ち、状態予測部２５までの処理では、各ファクタmの状態変化の認識及び予測を行うことはできるが、各ファクタmが住宅内の各場所に設置されたどの電気機器１２に対応するかは特定できない。機器特定部２６は、各ファクタと住宅内の電気機器１２とを対応付ける。例えば、機器特定部２６は、電気機器１２の製品の種類ごと（洗濯機、冷蔵庫、エアコン、など）に代表的な電流波形パターンを記憶し、その代表的な電流波形パターンと、実際に取得されたファクタの観測データYとを比較して、電気機器１２（の種類）を特定する。また例えば、機器特定部２６は、不図示の操作部において、各ファクタに対するユーザによる電気機器１２の種類の入力操作を受け付け、入力された電気機器１２の種類（製品）とファクタの観測データYとを対応付ける。
【００４７】
表示部２７は、液晶ディスプレイ等により構成され、パラメータ推定部２２による推定結果に基づいて、機器特定部２６により特定された電気機器１２の稼働状況を表示する。
【００４８】
以下、電気機器推定装置１を構成する各部の詳細について説明する。
【００４９】
[センサ部２１の詳細構成例]
図６は、センサ部２１の詳細構成例を示すブロック図である。
【００５０】
センサ部２１は、住宅内の各場所に設置された複数の電気機器１２の電流及び電圧の合計値を測定する電流測定部４１及び電圧測定部４２、電流波形切り出し部４３、並びに、時系列保持部４４により構成される。
【００５１】
図７は、単相三線式の配線により複数の電気機器１２に電気が供給される場合の、電流測定部４１と電圧測定部４２の例である。
【００５２】
電流測定部４１は、例えば、２つのクランプ型電流計５４及び５５により構成され、分電盤１１内の大元のアンペアブレーカの２次側に接続された２つの電圧線（ケーブル）５１及び５３にクランプされ、そこに流れる電流を計測する。
【００５３】
電圧測定部４２は、電圧計５６により構成され、所定の電線間の電圧を測定する。電圧は、電流の位相を電圧の位相に同期させるために測定される。従って、電圧測定部４２は、電圧線５１と中性線５２間の電圧、または、中性線５２と電圧線５３間の電圧のいずれか一方を計測すればよい。図７の例では、電圧計５６が電圧線５１と中性線５２間の電圧を測定している。
【００５４】
図６に戻り、電流波形切り出し部４３は、時々刻々と得られる電流値の時系列データ（電流波形）を時間方向に間引いて代表サンプルを生成し、時刻tの観測データYtとして時系列保持部４４に供給する。
【００５５】
図８は、電流測定部４１により測定された電流値の時系列データの例を示している。
【００５６】
図８上段の波形は、所定の期間内に、電流測定部４１により測定された電流値の時系列データであり、振幅の大きさが消費電流の大きさを表す。
【００５７】
図８中段の波形は、図８上段に示した時系列データのうち、４周期分の時系列データを拡大して示したものであり、さらに、図８下段の波形は、中段の４周期分の時系列データのうちの１周期分を拡大して示したものである。
【００５８】
電流波形切り出し部４３は、最初に、電流波形の位相を、電圧の位相に同期させる。具体的には、電流波形切り出し部４３は、電圧値のゼロ交差の瞬間を、電流波形の位相０として扱う。次に、電流波形切り出し部４３は、所定の時間間隔で、１周期分の電流波形をサンプリングする。全ての周期をサンプリングするとデータ量が膨大となってしまうからである。例えば、１秒間に１周期の電流波形を１サンプルとして取得するようにした場合、商用電源周波数５０Hzの地域においては、１／５０にデータ量を削減することができる。
【００５９】
商用電源周波数５０Hzの地域において、１回のサンプリングを５０kHzでサンプリングすることとすると、図８下段に示されるように、１回のサンプリングで１０００個の電流値が得られる。電流波形切り出し部４３は、１回のサンプリングで得られた１０００個の電流値を１０００次元の電流値ベクトルで表し、現在の時刻tにおける観測データYtとして、時系列保持部４４に供給する。従って、観測データYtは、D次元（D＝１０００）の実数ベクトルとなる。
【００６０】
図６に戻り、時系列保持部４４は、電流波形切り出し部４３から順次供給される観測データYtを記憶し、保持する。これにより、時系列保持部４４には、１サンプルが１０００次元の電流値ベクトルで表される観測データ｛Y１，Y２，Y３，・・・Yt，・・YT｝が格納される。
【００６１】
なお、上述した例では、単純に１秒間に１回サンプリングすることで、観測データを時間方向に間引いたが、観測データの間引き方はこれに限られない。例えば、５０サンプルの電流値ベクトルの平均ベクトルを計算し、その計算結果を１サンプルの観測データYtとしてもよい。また、５０サンプルの電流値ベクトルのうち、５０サンプルの平均ベクトルに最も近い電流値ベクトルを、１サンプルの観測データYtとしてもよい。あるいは、５０サンプルの電流値ベクトルをK-means法でクラスタリングし、最もサンプル数の多いクラスタの平均ベクトル、または、５０サンプルの電流値ベクトルのベクトル間距離が最小となる２サンプルの平均ベクトルを、１サンプルの観測データYtとしてもよい。
【００６２】
時系列保持部４４に格納されている観測データ｛Y１，Y２，Y３，・・・Yt，・・YT｝に基づいて、パラメータ推定部２２は、複数の電気機器１２の稼働状態をモデル化したFactorial HMMのモデルパラメータを推定する。
【００６３】
[Factorial HMMについて]
次に、パラメータ推定部２２の詳細について説明するが、初めに、Factorial HMMにおけるモデルパラメータの推定について説明する。
【００６４】
観測データ｛Y１，Y２，Y３，・・・，Yt，・・・，YT｝に対する隠れ状態を｛S１，S２，S３，・・・，St，・・・ST｝とすると、隠れ状態Sｔと観測データYｔの同時確率は、次式（１）で与えられる。
【数１】
0005598200.tif 000002
【００６５】
式（１）中の、P(S１)は初期確率、P(St｜St-1)は状態遷移確率、P(Yt｜St)は観測確率を表し、それぞれ、式（２）、式（３）、式（４）で計算される。
【数２】
0005598200.tif 000003
【００６６】
上述したように複数のファクタが１つの電気機器１２に対応する場合もあるが、最もシンプルな例として１つのファクタが１つの電気機器１２に対応するものとして、Factorial HMMにおけるモデルパラメータの推定について説明する。１つのファクタが１つの電気機器１２に対応するとした場合の、ファクタmに対応する電気機器１２を、m番目の電気機器１２とも称する。
【００６７】
式（２）乃至式（４）中のSt(m)は、時刻ｔにおけるm番目の電気機器１２の状態（オン、オフ、強運転、弱運転など）を表し、m番目の電気機器１２の状態数がKであるとすると、St(m)は、K次元の縦ベクトル（K行１列のベクトル）で構成される。
【００６８】
式（２）の初期確率P(S１)は、M個のπ(m)の掛け算で計算される。π(m)は、m番目の電気機器１２の初期状態確率を表し、K次元の縦ベクトルである。
【００６９】
式（３）の状態遷移確率P(St｜St-1)は、M個のA(m)の掛け算で計算される。A(m)は、例えば、オンからオフへの切り替わりやすさなどに対応する、m番目の電気機器１２の状態遷移確率を表し、K行K列（K×K）の正方行列で構成される。
【００７０】
式（４）の観測確率P(Yt｜St)は、観測平均μt，共分散行列Cの多変量正規分布で計算される。式（４）において、ダッシュ（’）は転置を表し、右上の”−１”は逆数を表す。また、｜C｜はCの絶対値を表す。
【００７１】
式（４）のW(m)は、m番目の電気機器１２が消費する電流波形のパターンに対応する観測確率P(Yt｜St)のパラメータである。電気機器１２の状態ごとに電流波形のパターンは異なるため、W(m)は、観測データの次元数Dを行数、状態数Kを列数としたD行K列（D×K）の行列となる。
【００７２】
μtは、時刻ｔにおける観測平均を表し、行列W(m)の状態St(m)に対応する列要素をM個足し合わせたものとなる。換言すれば、μtは、全ての電気機器１２の状態に応じた電流値を足し合わせたものに相当する。従って、観測平均μtが時刻ｔにおける観測データYtに近ければ、モデルパラメータが尤もらしいということになる。共分散行列Cは、電流パターンに乗るノイズの強度に対応し、全時刻、全ての電気機器１２で共通とされる。
【００７３】
以上より、パラメータ推定部２２によるFactorial HMMのモデルパラメータの推定とは、具体的には、m番目の電気機器１２の初期状態確率π(m)、状態遷移確率A(m)、観測確率のパラメータW(m)、及び共分散行列Cを求めることである。図９に、Factorial HMMのモデルパラメータφ＝｛π(m)，A(m)，W(m)，C｝とNILMにおける事象との対応関係を示す。
【００７４】
上述した文献Xには、モデルパラメータφ＝｛π(m)，A(m)，W(m)，C｝の推定方法として、４つの推定手法が開示されている。具体的には、１）厳密な推定 (exact inference)、２）ギブスサンプリングを使った近似推定 (inference using Gibbs sampling)、３）完全に分解した変分法による近似推定 (completely factorized variational inference)、４）構造を残した変分法による近似推定 (structured variational inference)が開示されている。
【００７５】
４つの推定手法は、いずれもEMアルゴリズムを用いた手法であり、パラメータ推定部２２が行う推定処理のベースとして採用することができる。「１）厳密な推定」は、演算量が多いが高精度で安定した推定が可能であるため、電気機器１２の数が少ない小規模のシステムや開発段階で利用することは有益である。実行速度、実行時のメモリ使用量、及び実行精度のそれぞれについて総合的に判断すると、「４）構造を残した変分法による近似推定」が最も有効である。
【００７６】
以下では、４つの推定手法のうち、「３）完全に分解した変分法による近似推定」と、「４）構造を残した変分法による近似推定」について、パラメータ推定部２２が行うとして説明する。
【００７７】
[完全に分解した変分法によるパラメータ推定処理]
図１０は、「３）完全に分解した変分法による近似推定」である、完全に分解した変分法によるパラメータ推定処理のフローチャートである。
【００７８】
初めに、ステップＳ１において、パラメータ推定部２２は、パラメータ推定処理における作業用変数などを初期化する初期化処理を行う。具体的には、パラメータ推定部２２は、全ての時刻t及びファクタm（t=1,・・・，T，m=1,・・・，M）について、変分パラメータθt(m)、ファクタmの観測確率のパラメータW(m)、共分散行列C、及び、状態遷移確率Ai,j(m)を初期化する。変分パラメータθt(m)と状態遷移確率Ai,j(m)には、初期値として1/Kが代入される。ファクタmの観測確率のパラメータW(m)には、初期値として所定の乱数が代入される。共分散行列Cの初期値には、C=aI（aは任意の実数、IはD行D列（D×D）の単位行列）が設定される。
【００７９】
ステップＳ２において、パラメータ推定部２２は、状態確率の推定を行うEステップ処理を実行する。Eステップ処理の詳細は、図１１を参照して後述する。
【００８０】
ステップＳ３において、パラメータ推定部２２は、遷移と観測のパラメータを推定するMステップ処理を実行する。Mステップ処理の詳細は、図１２を参照して後述する。
【００８１】
ステップＳ４において、パラメータ推定部２２は、モデルパラメータの収束条件を満たしたかを判定する。例えば、パラメータ推定部２２は、ステップＳ２乃至４の処理の繰り返し回数が予め設定した所定の回数に到達した場合、または、モデルパラメータの更新による状態尤度の変化量が所定値以内である場合に、モデルパラメータの収束条件を満たしたと判定する。
【００８２】
ステップＳ４で、モデルパラメータの収束条件をまだ満たしていないと判定された場合、処理はステップＳ２に戻り、ステップＳ２乃至Ｓ４の処理が繰り返される。
【００８３】
一方、ステップＳ４で、モデルパラメータの収束条件を満たしたと判定された場合、パラメータ推定部２２は、パラメータ推定処理を終了する。
【００８４】
[Eステップ処理の詳細フローチャート]
次に、図１１のフローチャートを参照して、図１０のステップＳ２として実行されるEステップ処理の詳細について説明する。
【００８５】
Eステップ処理では、初めに、ステップＳ１１において、パラメータ推定部２２は、ファクタに対応する変数mに１を代入する。
【００８６】
ステップＳ１２において、パラメータ推定部２２は、テンポラリ変数△(m)を、次式（５）により求める。
【数３】
0005598200.tif 000004
式（５）のdiag(・)は、（・）内の行列の対角成分をベクトルとして取り出す関数である。
【００８７】
ステップＳ１３において、パラメータ推定部２２は、式（６）と式（７）により、全ての時刻t（t=1,・・・，T）についてテンポラリ変数ηt(m)と変分パラメータθt(m)を求める。
【数４】
0005598200.tif 000005
式（６）のサメーション（Σ）は、ｌをm以外の１からMとしたときの和を意味する。また、式（７）のSoftmax関数は、式（８）で表される処理を行う関数である。
【数５】
0005598200.tif 000006
【００８８】
ステップＳ１４において、パラメータ推定部２２は、変数mが状態変数の数Mと等しいか、即ち、全てのファクタ１乃至Mについて変分パラメータθt(m)を求めたかを判定する。ステップＳ１４で、変数mが状態変数の数Mと等しくないと判定された場合、処理はステップＳ１５に進み、パラメータ推定部２２は、変数mを１だけ増加させた後、処理をステップＳ１２に戻す。これにより、更新後の変数mについて、ステップＳ１２乃至Ｓ１４の処理が繰り返される。
【００８９】
一方、ステップＳ１４で、変数mが状態変数の数Mと等しいと判定された場合、即ち、全てのファクタ１乃至Mについて変分パラメータθt(m)が得られた場合、処理はステップＳ１６に進み、パラメータ推定部２２は、変分パラメータθt(m)の収束条件を満たしたかを判定する。ステップＳ１６では、例えば、繰り返し回数が予め設定した所定の回数に到達した場合に、変分パラメータθt(m)の収束条件を満たしたと判定される。
【００９０】
ステップＳ１６で、変分パラメータθt(m)の収束条件を満たしていないと判定された場合、処理はステップＳ１１に戻る。そしてステップＳ１１乃至Ｓ１６が再度実行されることにより、再度、変分パラメータθt(m)が計算（更新）される。
【００９１】
一方、ステップＳ１６で、変分パラメータθt(m)の収束条件を満たしたと判定された場合、処理はステップＳ１７に進み、パラメータ推定部２２は、再度、ファクタに対応する変数mに１を代入する。
【００９２】
そして、ステップＳ１８において、パラメータ推定部２２は、期待値変数<St(m)>，<S’t(m)>，<St(m)S’t(m)>を、次式（９）乃至式（１１）により求める。
【数６】
0005598200.tif 000007
【００９３】
さらに、ステップＳ１９において、パラメータ推定部２２は、変数ｎを１乃至Mとした期待値変数<St(m)S’t(n)>を、次式（１２）により求める。
【数７】
0005598200.tif 000008
【００９４】
ステップＳ２０において、パラメータ推定部２２は、変数mが状態変数の数Mと等しいか、即ち、全てのファクタ１乃至Mについて、期待値変数<St(m)>，<S’t(m)>，<St(m)S’t(m)>、及び<St(m)S’t(m)>を求めたかを判定する。
【００９５】
ステップＳ２０で、変数mが状態変数の数Mと等しくないと判定された場合、処理はステップＳ２１に進み、パラメータ推定部２２は、変数mを１だけ増加させた後、処理をステップＳ１８に戻す。これにより、更新後の変数mについて、ステップＳ１８乃至Ｓ２０の処理が再度繰り返される。
【００９６】
一方、ステップＳ２０で、変数mが状態変数の数Mに等しいと判定された場合、即ち、全てのファクタ１乃至Mについて期待値変数<St(m)>，<S’t(m)>，<St(m)S’t(m)>、及び<St(m)S’t(m)>が得られた場合、Eステップ処理は終了する。
【００９７】
[Mステップ処理の詳細フローチャート]
Eステップ処理が終了すると、図１０のステップＳ３におけるMステップ処理が実行される。図１２は、ステップＳ３におけるMステップ処理の詳細を説明するフローチャートである。
【００９８】
Mステップ処理では、ステップＳ３１において、パラメータ推定部２２は、全てのファクタm＝１乃至Mについて、初期状態確率π(m)を、次式（１３）により求める。
【数８】
0005598200.tif 000009
【００９９】
ステップＳ３２において、パラメータ推定部２２は、全てのファクタmについて、状態Sj(m)から状態Si(m)への状態遷移確率Ai,j(m)を、次式（１４）により求める。
【数９】
0005598200.tif 000010
ここで、St-1,j(m)は、遷移前の状態Sj(m)が時刻ｔ−１の状態変数St-1(m)であり、St,i(m)は、遷移後の状態Si(m)が時刻ｔの状態変数St(m)であることを表す。
【０１００】
ステップＳ３３において、パラメータ推定部２２は、観測確率のパラメータWを、次式（１５）により求める。
【数１０】
0005598200.tif 000011
【０１０１】
式（１５）において、観測確率のパラメータWは、D行K列（D×K）のパラメータW(m)を、列方向に全てのファクタmについてM個連結した、D行MK列（D×MK，MKはMとKの積）の行列を表す。従って、ファクタmの観測確率のパラメータW(m)は、観測確率のパラメータWを列方向に分解することで得られる。また、式（１５）におけるpinv(・)は、疑似逆行列を求める関数である。
【０１０２】
ステップＳ３４において、パラメータ推定部２２は、共分散行列Cを、次式（１６）により求める。
【数１１】
0005598200.tif 000012
【０１０３】
上述したステップＳ３１乃至Ｓ３４により、Factorial HMMのモデルパラメータφが求められ（更新され）、Mステップ処理が終了して、図１０のステップＳ４の処理が進められる。
【０１０４】
[構造を残した変分法によるパラメータ推定処理]
次に、「４）構造を残した変分法による近似推定」である、構造を残した変分法によるパラメータ推定処理について説明する。図１３は、構造を残した変分法によるパラメータ推定処理のフローチャートである。
【０１０５】
初めに、ステップＳ４１において、パラメータ推定部２２は、パラメータ推定処理における作業用変数などを初期化する初期化処理を行う。具体的には、パラメータ推定部２２は、全ての時刻t及びファクタm（t=1,・・・，T，m=1,・・・，M）について、変分パラメータｈt(m)、期待値変数<St(m)>、ファクタmの観測確率のパラメータW(m)、共分散行列C、及び、状態遷移確率Ai,j(m)を初期化する。変分パラメータｈt(m)には初期値として1/Kが代入され、それ以外については、上述したステップＳ１の初期化処理と同様の値が代入される。
【０１０６】
続くステップＳ４２乃至Ｓ４４の処理は、図１０のステップＳ２乃至Ｓ４と、それぞれ同様であるので、その説明は省略する。
【０１０７】
[Eステップ処理の詳細フローチャート]
次に、図１４のフローチャートを参照して、図１３のステップＳ４２として実行されるEステップ処理の詳細について説明する。
【０１０８】
構造を残した変分法によるパラメータ推定処理のEステップ処理では、初めに、ステップＳ６１において、パラメータ推定部２２は、ファクタに対応する変数mに１を代入する。
【０１０９】
そして、ステップＳ６２において、パラメータ推定部２２は、テンポラリ変数△(m)を、上述した式（５）により求める。即ち、ステップＳ６１及びＳ６２の処理は、上述した図１１のステップＳ１１及びＳ１２と同様である。
【０１１０】
ステップＳ６３において、パラメータ推定部２２は、式（１７）と式（１８）により、全ての時刻t（t=1,・・・，T）についてテンポラリ変数ηt(m)と変分パラメータｈt(m)を求める。
【数１２】
0005598200.tif 000013
【０１１１】
ステップＳ６４において、パラメータ推定部２２は、変数mが状態変数の数Mと等しいか、即ち、全てのファクタ１乃至Mについて変分パラメータｈt(m)を求めたかを判定する。ステップＳ６４で、変数mが状態変数の数Mと等しくないと判定された場合、処理はステップＳ６５に進み、パラメータ推定部２２は、変数mを１だけ増加させた後、処理をステップＳ６２に戻す。これにより、更新後の変数mについて、ステップＳ６２乃至Ｓ６４の処理が繰り返される。
【０１１２】
一方、ステップＳ６４で、変数mが状態変数の数Mと等しいと判定された場合、即ち、全てのファクタ１乃至Mについて変分パラメータｈt(m)が得られた場合、処理はステップＳ６６に進み、パラメータ推定部２２は、変分パラメータｈt(m)の収束条件を満たしたかを判定する。ステップＳ６６では、例えば、繰り返し回数が予め設定した所定の回数に到達した場合、パラメータ推定部２２は、変分パラメータｈt(m)の収束条件を満たしたと判定する。
【０１１３】
ステップＳ６６で、変分パラメータｈt(m)の収束条件を満たしていないと判定された場合、処理はステップＳ６１に戻る。そしてステップＳ６１乃至Ｓ６６が再度実行されることにより、再度、変分パラメータｈt(m)が計算（更新）される。
【０１１４】
一方、ステップＳ６６で、変分パラメータｈt(m)の収束条件を満たしたと判定された場合、処理はステップＳ６７に進み、パラメータ推定部２２は、変分パラメータｈt(m)と状態遷移確率Ai,j(m)を用いたフォーワードバックワードアルゴリズムを用いて、期待値変数<St(m)>及び<St-1(m)S’t(m)>を求める。
【０１１５】
ステップＳ６８において、パラメータ推定部２２は、変数ｎを１乃至Mとした期待値変数<St(m)S’t(n)>を次式（１９）により求め、Eステップ処理を終了して、図１３に戻る。
【数１３】
0005598200.tif 000014
【０１１６】
構造を残した変分法によるパラメータ推定処理において、Eステップ処理の次に実行されるMステップ処理は、図１２を参照して説明した、完全に分解した変分法によるパラメータ推定処理のMステップ処理と同様であるので、その説明は省略する。
【０１１７】
[パラメータ推定部２２によるFactorial HMMの手法]
以上、文献Xに開示されているFactorial HMMのモデルパラメータφを求める４つの推定手法のうち、完全に分解した変分法によるパラメータ推定処理と、構造を残した変分法によるパラメータ推定処理について説明した。
【０１１８】
Factorial HMM自体は、汎用の時系列データのモデル化手法であり、NILM以外の様々な問題に適用可能なものである。しかしながら、Factorial HMMは汎用モデルであるが故に、Factorial HMMをNILMにそのまま適用すると、次の２つの問題が生じてしまう。
【０１１９】
１）物理的に存在し得ない「負の消費電力」を持つファクタが発生する。これは、Factorial HMMのパラメータの自由度が大きすぎるため、Factorial HMMのパラメータW(m)が負の値をも取り得ることに対応する。
２）各ファクタと住宅内に設置された各電気機器１２とが１対１に対応しないことがある。換言すれば、１つの電気機器１２が複数のファクタと対応する場合がある。これは、Factorial HMMの表現力が高いために、一つの観測データを二通り以上の方法で説明できてしまうからである。
【０１２０】
そこで、本発明を適用した電気機器推定装置１のパラメータ推定部２２では、NILMに適用するため、ベースとなる文献Xに開示されているFactorial HMMに対して、以下に示す改良を施した手法が採用されている。
【０１２１】
パラメータ推定部２２は、Factorial HMMの改良として、第１に、Factorial HMMにおける観測確率のパラメータW(m)が非負行列であるという制約（非負制約）を付加する。供給電源は交流であるため、観測データYとなる電流値は正値と負値を交互に取る。そこで、パラメータ推定部２２は、次のいずれかの方法を前処理として行うことにより、観測データYとなる電流値を正のみの値に変換する。具体的には、パラメータ推定部２２は、１）電流値の絶対値を観測データYとして用いる、２）正値を取る半周期分のみを観測データYとして用いる、３）電力値（電流値と電圧値の積）を観測データYとして用いる。
【０１２２】
さらに、パラメータ推定部２２は、上述した式（１５）を用いず、次式（２０）で表される目的関数を最小化する観測確率のパラメータWを、制約付き２次計画法により求める。
【数１４】
0005598200.tif 000015
【０１２３】
即ち、パラメータ推定部２２は、上述した式（１５）に代えて、制約条件Wvertical≧０の下で、目的関数としての、ＦとWverticalの積ベクトルとｇベクトルとの差の絶対値の二乗を最小化するWverticalを求める。式（２０）は、Factorial HMMの全てのファクタについての、パラメータW(m)が表す電流波形のパターンと、観測データYが表す前記消費電流の合計値のパターンとの誤差を最小とするWverticalを求めることを意味する。なお、reshape_to_vertical_vector(・)は、行列の列成分を縦に連結して、例えば、ａ行ｂ列（ａ×ｂ）の行列を、ａｂ行１列（ab×１，abはａとｂの積）のベクトルに変形する操作を表す関数である。
【０１２４】
求めるべきパラメータWverticalには、
Wvertical＝reshape_to_vertical_vector(W)
という関係があるため、式（２０）で求められたａｂ行１列のWverticalを、ａ行ｂ列に変換することで、観測確率のパラメータWが得られる。
【０１２５】
式（２０）の最適化処理の実行には、文献「Lawson, C. L. and R. J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974, Chapter 23.」を基に実装したFortranのコードnnls.fや、そのC言語版であるnnls.cが広く知られている。また、Matlab（登録商標）やPython （登録商標）などでも、式（２０）の最適化処理を簡単に実行するツールが提供されている。
【０１２６】
パラメータ推定部２２は、Factorial HMMの改良として、第２に、状態数Kを、オン、オフのみに対応する２（K＝２）に固定し、観測確率のパラメータWの第１列目または第２列目のいずれか一方を常に０ベクトルに設定する。これは、電気機器１２の状態がオンとオフの２状態しかないこと、及び、状態がオフのときは、消費電力が常に０であること、を仮定したものである。
【０１２７】
このオンオフ２状態の制約条件の利点は、Factorial HMMのパラメータ推定結果が現実の電気機器１２の動作モードに対応するようになり、Factorial HMMのパラメータの学習の安定性が向上することである。白熱灯などの、いくつかのシンプルな電気機器１２は、このオンオフ２状態の制約条件を完全に満たすが、動作モードが複数ある電気機器１２であっても、現実の電気機器１２の動作モードに、より対応するようになる。この制約条件を付さないと、Factorial HMMのパラメータ推定結果が現実の電気機器の動作モードと対応しない一つの理由としては、Factorial HMMは非常に表現力豊かな手法であるため、現実には対応しないが高い尤度を持つ局所解がいくつも存在してしまうことが考えられる。
【０１２８】
[パラメータ推定部２２によるMステップ処理]
図１５を参照して、パラメータ推定部２２によるMステップ処理について説明する。
【０１２９】
なお、Mステップ処理以外の処理については、図１０及び図１１、図１３及び図１４等を参照して説明した処理と同様である。但し、時系列保持部４４に保存されている観測データYが正のみの値に変換されたものである点、及び、パラメータW(m)の状態数Kが２（K＝２）で、第１列目または第２列目のいずれか一方が常に０ベクトルである点は異なる。
【０１３０】
Mステップ処理では、ステップＳ８１において、パラメータ推定部２２は、全てのファクタm＝１乃至Mについて、初期状態確率π(m)を、上述した式（１３）により求める。
【０１３１】
ステップＳ８２において、パラメータ推定部２２は、全てのファクタmについて、状態Sj(m)から状態Si(m)への状態遷移確率Ai,j(m)を、上述した式（１４）により求める。
【０１３２】
ステップＳ８３において、パラメータ推定部２２は、観測確率のパラメータWを、上述した式（１５）を用いず、式（２０）の制約付き２次計画法によって求める。
【０１３３】
ステップＳ８４において、パラメータ推定部２２は、共分散行列Cを、上述した式（１６）により求める。
【０１３４】
以上の処理により、パラメータ推定部２２によるMステップ処理が実行される。
【０１３５】
[同一機器判定部２４の処理]
次に、同一機器判定部２４の処理について説明する。
【０１３６】
図１６は、本来の状態数Kが３（K＝３）である電気機器１２、例えば、停止モード、弱運転モード、及び強運転モードの３つの動作モードを持つ扇風機を、オンオフ２状態の制約条件下で推定したときのパラメータ推定結果を示している。
【０１３７】
本来の状態数Kが３（K＝３）である電気機器１２をオンオフ２状態の制約条件下で推定すると、パラメータ推定結果は、図１６の可能性Aか、または、可能性Bに収束する。
【０１３８】
図１６では、オンオフ２状態として推定したときの推定結果としての状態１と状態２を分かりやすくするため、”ＯＮ”と”ＯＦＦ”と記述している。実際には、例えば、状態１がOFF（オフ）で状態２がON（オン）として、ONのときはSt(m)=[状態１，状態２]＝[0.0，1.0]、オフのときはSt(m)= [状態１，状態２]＝[0.9, 0.1]などのようになっている。
【０１３９】
可能性Aは、強運転モードを２つのファクタの同時生起として表現し、可能性Bは弱運転モードと強運転モードを２つのファクタに別々に表現している。
【０１４０】
可能性Aと可能性Bはどちらも、扇風機の稼動状態を正しくモデル化しているが、これら２つのファクタが２つの別々の電気機器１２を表現するのか、１つの電気機器１２の異なるモードを表現しているのかは、このままではわからない。しかし、１つの電気機器１２である場合は、以下の条件が常に成立することがわかっている。
【０１４１】
可能性Aの場合：
必要条件  ：  S(2)=ONのとき、必ずS(1)=ONになっている。
ヒント条件：  S(1)=ONからS(1)=OFFに遷移するとき、同時刻にS(2)=ONからS(2)=OFFに遷移することがある。
ヒント条件：  S(1)=OFFからS(1)=ONに遷移するとき、同時刻にS(2)=OFFからS(2)=ONに遷移することがある。
可能性Bの場合：
必要条件：  S(1)=ONのとき、必ずS(2)=OFFになっている。
必要条件：  S(2)=ONのとき、必ずS(1)=OFFになっている。
ヒント条件：  S(1)=ONからS(1)=OFFに遷移するとき、同時刻にS(2)=OFFからS(2)=ONに遷移することがある。
ヒント条件：  S(1)=OFFからS(1)=ONに遷移するとき、同時刻にS(2)=ONからS(2)=OFFに遷移することがある。
【０１４２】
従って、この条件を評価することで、状態S(1)と状態S(2)が同一の電気機器１２に起因していることを推定することが可能である。
【０１４３】
図１７は、同一機器判定部２４による同一機器判定処理のフローチャートを示している。
【０１４４】
初めに、ステップＳ１０１において、同一機器判定部２４は、Factorial HMMの各ファクタの全組み合わせを、同一機器候補として作成する。
【０１４５】
ステップＳ１０２において、同一機器判定部２４は、必要条件による同一機器候補の選別を行う。即ち、同一機器判定部２４は、各ファクタの全組み合わせのうち、必要条件を常に満たす組み合わせを、同一機器候補として選別する。
【０１４６】
ステップＳ１０３において、同一機器判定部２４は、ヒント条件による同一機器候補の選別を行う。即ち、同一機器判定部２４は、ステップＳ１０２で選別された組み合わせのうち、ヒント条件をZ回（Zは２以上の整数）以上満たしたことがある組み合わせを同一機器と判定する。同一機器判定部２４は、同一機器と判定したファクタについての情報をデータベース２３に保存して、処理を終了する。
【０１４７】
[電気機器推定装置１の全体処理]
図１８のフローチャートを参照して、電気機器推定装置１の処理全体について説明する。
【０１４８】
初めに、ステップＳ１２１において、センサ部２１は、観測データYとしての、住宅内の消費電流の合計値を測定し、取得された合計値に対して、正のみの値に変換する前処理を施してから、時系列保持部４４に保存する。
【０１４９】
ステップＳ１２２において、パラメータ推定部２２は、観測データ｛Y１，Y２，Y３，・・・Yt，・・YT｝に基づいて、Factorial HMMのモデルパラメータを推定するパラメータ推定処理を行う。そして、推定結果としてのモデルパラメータ（の値）が、データベース２３に保存される。
【０１５０】
ステップＳ１２２のパラメータ推定処理としては、例えば、図１０の完全に分解した変分法によるパラメータ推定処理や図１３の構造を残した変分法によるパラメータ推定処理が行われる。ただし、図１０のステップＳ３及び図１３のステップＳ４３のMステップ処理としては、図１５のMステップ処理が実行される。
【０１５１】
ステップＳ１２３において、同一機器判定部２４は、図１７を参照して説明した同一機器判定処理を行う。
【０１５２】
ステップＳ１２４において、状態予測部２５は、現在時刻から所定時間経過後の未来のファクタmの状態を予測する予測処理を行う。具体的には、状態予測部２５は、ファクタmが現在時刻Tから所定時間L経過後の未来の時刻T+Lにどの状態になっているかを表す、時刻T+Lにおけるファクタmの状態確率ST+L(m)を、次式（２１）により求める。
【数１５】
0005598200.tif 000016
【０１５３】
状態予測部２５は、予測結果を、例えば、図示せぬ関連機器制御部等に出力して、処理を終了する。
【０１５４】
以上、説明した電気機器推定装置１の電気機器推定処理によれば、住宅内の各場所に設置された複数の電気機器１２の消費電流の合計値でなる観測データYをFactorial HMMによりモデル化し、モデルパラメータを求める。このFactorial HMMには、NILMに適用するため、従来のFactorial HMMに対して改良が施されている。具体的には、電気機器推定装置１は、観測確率のパラメータW(m)を非負の値となるように変換し、制約付き２次計画法により求める。また、電気機器推定装置１は、状態変数S(１)乃至S(M)の状態数Kを、オン、オフのみに対応する２（K＝２）に固定し、観測確率のパラメータWの第１列目または第２列目のいずれか一方を常に０ベクトルとする。換言すれば、Factorial HMMの各状態が取り得る状態数を２とし、Factorial HMMのファクタmの電流波形のパターンに対応する観測確率のパラメータW(m)が非負であって、状態数に対応する第１列目または第２列目のいずれか一方が常に０であるという制約条件の下で、Factorial HMMが、時系列データが表す消費電流の合計値のパターンを説明する度合いである尤度関数を最大化することにより、モデルパラメータとしての観測確率のパラメータW(m)が求められる。これにより、事前のデータベースが不要で、Factorial HMMのモデルパラメータを容易に求めることができる。即ち、パラメータの算出及び調整が簡単であり、かつ、事前のデータベースが不要な手法を確立することができる。
【０１５５】
上述した一連の処理は、ハードウエアにより実行することもできるし、ソフトウエアにより実行することもできる。一連の処理をソフトウエアにより実行する場合には、そのソフトウエアを構成するプログラムが、コンピュータにインストールされる。ここで、コンピュータには、専用のハードウエアに組み込まれているコンピュータや、各種のプログラムをインストールすることで、各種の機能を実行することが可能な、例えば汎用のパーソナルコンピュータなどが含まれる。
【０１５６】
図１９は、上述した一連の処理をプログラムにより実行するコンピュータのハードウエアの構成例を示すブロック図である。
【０１５７】
コンピュータにおいて、CPU（Central Processing Unit）１０１，ROM（Read Only Memory）１０２，RAM（Random Access Memory）１０３は、バス１０４により相互に接続されている。
【０１５８】
バス１０４には、さらに、入出力インタフェース１０５が接続されている。入出力インタフェース１０５には、入力部１０６、出力部１０７、記憶部１０８、通信部１０９、及びドライブ１１０が接続されている。
【０１５９】
入力部１０６は、キーボード、マウス、マイクロホンなどよりなる。出力部１０７は、ディスプレイ、スピーカなどよりなる。記憶部１０８は、ハードディスクや不揮発性のメモリなどよりなる。通信部１０９は、ネットワークインタフェースなどよりなる。ドライブ１１０は、磁気ディスク、光ディスク、光磁気ディスク、或いは半導体メモリなどのリムーバブル記録媒体１１１を駆動する。
【０１６０】
以上のように構成されるコンピュータでは、CPU１０１が、例えば、記憶部１０８に記憶されているプログラムを、入出力インタフェース１０５及びバス１０４を介して、RAM１０３にロードして実行することにより、上述した一連の処理が行われる。
【０１６１】
コンピュータ（CPU１０１）が実行するプログラムは、例えば、パッケージメディア等としてのリムーバブル記録媒体１１１に記録して提供することができる。また、プログラムは、ローカルエリアネットワーク、インターネット、デジタル衛星放送といった、有線または無線の伝送媒体を介して提供することができる。
【０１６２】
コンピュータでは、プログラムは、リムーバブル記録媒体１１１をドライブ１１０に装着することにより、入出力インタフェース１０５を介して、記憶部１０８にインストールすることができる。また、プログラムは、有線または無線の伝送媒体を介して、通信部１０９で受信し、記憶部１０８にインストールすることができる。その他、プログラムは、ROM１０２や記憶部１０８に、あらかじめインストールしておくことができる。
【０１６３】
なお、本明細書において、フローチャートに記述されたステップは、記載された順序に沿って時系列的に行われる場合はもちろん、必ずしも時系列的に処理されなくとも、並列に、あるいは呼び出しが行われたとき等の必要なタイミングで実行されてもよい。
【０１６４】
なお、本明細書において、システムとは、複数の装置により構成される装置全体を表すものである。
【０１６５】
本発明の実施の形態は、上述した実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において種々の変更が可能である。
【符号の説明】
【０１６６】
１  電気機器推定装置，  １１  分電盤，  １２  電気機器，  ２１  センサ部，  ２２  パラメータ推定部，  ２３  データベース，  ２４  同一機器判定部，  ２５  状態予測部，  ４１  電流測定部，  ４２  電圧測定部，  ４３  電流波形切り出し部，  ４４  時系列保持部
【図１】
0005598200.tif 000017
【図２】
0005598200.tif 000018
【図３】
0005598200.tif 000019
【図４】
0005598200.tif 000020
【図５】
0005598200.tif 000021
【図６】
0005598200.tif 000022
【図７】
0005598200.tif 000023
【図８】
0005598200.tif 000024
【図９】
0005598200.tif 000025
【図１０】
0005598200.tif 000026
【図１１】
0005598200.tif 000027
【図１２】
0005598200.tif 000028
【図１３】
0005598200.tif 000029
【図１４】
0005598200.tif 000030
【図１５】
0005598200.tif 000031
【図１６】
0005598200.tif 000032
【図１７】
0005598200.tif 000033
【図１８】
0005598200.tif 000034
【図１９】
0005598200.tif 000035