CN106170708A - 个体电设备工作状态估计装置及其方法 - Google Patents

Abstract

个体电设备工作状态估计装置具有：总使用电力测定单元(5)；时间变动计算单元(6a)，其计算总使用电力的时间性变动(跳跃电力)；电设备机型确定单元(6b)，其根据跳跃电力来判定是否要作为新电设备追加到现有的机型列表、或删除多余的机型；对应似然度估计单元(6c)，其估计关于升和降的跳跃与哪个电设备的启动/关闭的事件的发生相对应的似然度；以及个体电设备工作状态估计单元(6d)，其估计在输入了似然度时的个体电设备的工作状态的变化，对当前的工作状态的估计进行更新，根据使用工作状态的估计结果和标准电力的估计结果而得到的总使用电力的预测值以及总使用电力的测定数据来优化个体电设备的机型、标准电力以及工作状态的各估计，动态地估计个体电设备的工作概率。

Description

个体电设备工作状态估计装置及其方法

技术领域

本发明涉及一种能够将多个电设备中的正在工作的电设备的种类及其工作状态同时实时地进行监视的个体电设备工作状态估计装置及其方法。

背景技术

为了实现电力的高效利用，期望的是，针对电设备能够个别地掌握各自的使用状况。在该情况下，如果对各电设备安装传感器，则能够个别地测定各个电设备的使用状况(每个电设备的工作/非工作的状态、其电力使用量)，基于该测定结果来掌握电设备的使用状况，由此能够有助于节能化、电力供给控制等。

但是，若像这样对各电设备分别设置传感器，则不仅传感器数量增加，而且需要这些传感器的安装成本、用于构建收集来自传感器的测定信号的收集系统(无线等)的成本，从而无法避免高成本化。

因此，开发了一种即使不对各电设备分别设置传感器也能够监视各电设备的使用状况的技术，以作为廉价的装置在一般家庭等中也能够利用。

关于这种以往的个体电设备工作状态估计装置及其方法，已知专利文献1所记载的技术。

在该专利文献1所记载的个体电设备工作状态估计装置中，进行非侵入式电力使用监视(NIALM：Non-intrusive Appliance Load Monitoring)，即基于正在使用的电设备的总使用电力的电特性，利用脉冲型波形的模式识别，来个别地测定电设备的电力使用、操作特性。

即，更具体地说，该以往的个体电设备工作状态估计装置具备：传感器单元(电压传感器、电流传感器)，其在设备外的位置与电路连接，测定设备的上述电路中的全负荷的电特性(包括电导、电纳的导纳)，生成与它们成正比的模拟信号；AD转换器，其根据测定出的电压、电流来运算电导、电纳，将上述模拟信号转换为数字信号；以及信号处理器，其根据数字信号将多个单元各自的全负荷分为它们各自的部分。此外，能够利用该处理结果来将电设备的工作状况个别地显示在显示器上。

专利文献1：美国专利第4,858,141号公报

发明内容

发明要解决的问题

然而，在上述以往的个体电设备工作状态估计装置中存在如下面说明的那样的问题。

在上述以往技术中，在其应用时，需要事先收集哪个家庭持有几台具有什么样的标准电力的什么样的电设备之类的信息，或者在每当存在废弃的电设备、新追加的电设备的情况下都需要更新其信息。因此，耗费很大工夫、成本，而且无法指望在一般家庭中由外行人输入或更新这些信息，从而存在现实中非常难以应用之类的问题。

本发明是着眼于上述问题而完成的，其目的在于提供一种即使不对各个电设备设置传感器、而且事先不具有存在几台什么种类的电设备之类的信息也能够实时地对电设备的种类的确定以及它们的工作状态进行监视的个体电设备工作状态估计装置及其方法。

用于解决问题的方案

为了达到该目的，在本发明的个体电设备工作状态估计装置及其方法中，通过统计性逆估计的方法来构建根据总使用电力实时地估计个体电设备的机型及其工作状态的动态算法。因此，执行以下算法：根据总使用电力的时间性变动，对发生频度高的变动(跳跃电力)进行聚类来估计每个个体电设备(机型)的标准电力，根据机型的簇和测量出的跳跃来实时地估计电设备的确定及其工作状况，并且根据需要来进行簇的变更。

即，本发明的个体电设备工作状态估计装置的特征在于，具备：

总使用电力测定单元，其测定多个电设备的总使用电力；

时间变动计算单元，其基于由总使用电力测定单元测定出的总使用电力来计算跳跃电力，该跳跃电力是该总使用电力的时间性变动；

电设备机型确定单元，其进行聚类，即，基于由时间变动计算单元计算出的跳跃电力来判定是否要作为新的电设备追加到现有的机型列表，在要作为新的电设备追加到现有的机型列表的情况下，按个体电设备重新计算标准电力，根据需要而从现有的机型列表删除多余的机型，在不追加新的电设备的情况下维持现有的机型列表；

对应似然度估计单元，其基于由时间变动计算单元计算出的跳跃电力来估计关于升和降的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度；以及

个体电设备工作状态估计单元，其估计在输入了由对应似然度估计单元得到的、在当前的个体电设备的工作状态下根据总使用电力的跳跃而得到的、与个体电设备的启动和关闭的事件发生有关的似然度时的个体电设备的工作状态的变化，对当前的个体电设备的工作状态的估计进行更新，根据使用该更新后的个体电设备的工作状态的估计结果和标准电力的估计结果所预测出的总使用电力的预测值以及总使用电力的测定数据，来优化个体电设备的机型的估计、标准电力的估计以及工作状态的估计，以测定出的总使用电力的时间变动的时间间隔来动态地估计时刻变化的个体电设备的工作概率。

另一方面，本发明的个体电设备工作状态估计方法的特征在于，具有以下步骤：

测定多个电设备的总使用电力；

基于所测定出的总使用电力来计算跳跃电力，该跳跃电力是该总使用电力的时间性变动；

进行聚类，即，基于所计算出的跳跃电力来判定是否要作为新的电设备追加到现有的机型列表，在要作为新的电设备追加到现有的机型列表的情况下，按个体电设备重新计算标准电力，根据需要而从现有的机型列表删除多余的机型，在不追加新的电设备的情况下维持现有的机型列表；

基于所计算出的跳跃电力来估计关于升和降的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度；以及

估计在输入了在当前的个体电设备的工作状态下根据总使用电力的跳跃而得到的、与个体电设备的启动和关闭的事件发生有关的似然度时的该个体电设备的工作状态的变化，对当前的个体电设备的工作状态的估计进行更新，根据使用该更新后的个体电设备的工作状态的估计结果和标准电力的估计结果所预测出的总使用电力的预测值以及总使用电力的测定数据，来优化个体电设备的机型的估计、标准电力的估计以及工作状态的估计，以测定出的总使用电力的时间变动的时间间隔来动态地估计时刻变化的个体电设备的工作概率。

发明的效果

在本发明的个体电设备工作状态估计装置及其方法中，即使不对各个电设备设置传感器、而且事先不具有存在几台什么种类的电设备之类的信息也能够实时地监视作为监视对象的各个电设备的使用状况。

附图说明

图1是表示在具备多个电设备的房屋内设置有本发明的实施例1的个体电设备工作状态估计装置的状态的图。

图2是表示实施例1的个体电设备工作状态估计装置的结构的功能框图。

图3是表示示出了在实施例1的个体电设备工作状态估计装置中执行的个体电设备工作状态估计的流程的考虑方法的流程图的图。

图4是用于说明在实施例1的个体电设备工作状态估计装置中执行的k-均值算法(k＝2的情况)的示意图。

图5是表示对在实施例1的个体电设备工作状态估计装置中使用的马尔可夫转换模型中跳跃为正的情况下的状态变化进行总结而得到的表的图。

图6是表示对在实施例1的个体电设备工作状态估计装置中执行的马尔可夫转换模型中跳跃为负的情况下的状态变化进行总结而得到的表的图。

图7是表示将在实施例1的个体电设备工作状态估计装置中执行的个体电设备工作状态估计的流程与数式一起表示的流程图的图。

图8是表示对实施例1的个体电设备工作状态估计装置的效果进行确认的实验中的作为实际的家庭环境数据的一部分的天气的数据的图。

图9是表示上述实验中的按总使用电力的1分钟的跳跃量的规模的频数分布的图。

图10是表示上述实验中的按总使用电力规模的在总电力使用量中的占有率的图。

图11是表示上述实验中的8天内的实际环境数据的实时处理算法的结果(从开始使用起11520分钟为止)的图。

图12是放大了图11的结果中的从开始使用到150分钟为止的部分的图。

图13是放大了图11的结果中的从150分钟到350分钟为止的部分的图。

图14是表示使用以手动开启和关闭个体电设备的模拟家庭环境下的4小时内的以1分钟为间隔的数据并应用本发明来估计个体电设备的种类和标准电力、并估计总使用电力和个体电设备的工作状态而得到的结果的图。

图15是表示估计个体电设备的种类和标准电力、并在考虑了个体电设备的种类的数量发生变化时的状态变量的传递的基础上估计总使用电力和个体电设备的工作状态而得到的结果的图。

具体实施方式

首先，本发明的个体电设备工作状态估计装置及其方法在下述方面与上述专利文献1所记载的技术不同。即，它们在以下方面不同：

(a)是否仅着眼于作为总使用电力的时间变动的跳跃(jump)、其中特别是各电设备的开启-关闭；

(b)是否采用k-均值法(k-Means Clustering)的算法等统计性方法；

(c)需要多大程度的关于什么样的电设备正在工作等的先验信息；

(d)是否建立了能够应对以1分钟为间隔的测量等、实际环境下的总使用电力的测量时间序列数据的动态模型化。

作为具备所有这些性质的算法，在本发明中使用逆估计算法，但是目前不存在满足上述全部性质的NIALM的解法算法。

下面，基于附图所示的实施例来详细说明本发明的实施方式。

实施例1

在实施例1中，说明将本实施例的个体电设备工作估计装置应用于一般房屋(家庭)的例子。

图1示意性地表示将如上所述的使用逆估计算法的本实施例1的个体电设备工作状态估计装置设置在具备作为该个体电设备工作状态估计装置的监视对象的多个电设备的房屋中的状态下的整体结构。此外，逆估计算法是指基于不完整的数据以统计方式复原出完整数据的方法。

在该图中，在房屋1内，多条布线2从配电盘8延伸并遍布在其周围，在这些布线2的适当位置处分别经由接头、插座等而连接或能够连接多个机型的电设备3a～3n(n为2以上的正整数)。作为这种电设备3a～3n，例如包括如冰箱、电视、吸顶灯、空调(经由配电盘8并通过相对于前者等独立的系统来连接)等那样能够开启/关闭的家庭用电设备，除此以外还包括如卫洗丽(Washlet，商品名)那样的不定期地产生负荷的电设备、如对讲机等那样始终开启的电设备。

在房屋1内的配电盘8与屋外的供给电线4之间设置有电力计5，对该房屋1中正在使用的总使用电力进行测定。另外，也有时通过配电盘8来测定每个断路器的电力。所述电力计5、配电盘8相当于本发明的总使用电力测定单元。下面，说明个体电设备工作状态估计装置利用由电力计5测定出的总使用电力的情况，但是也可以利用配电盘8来测定每个该配电盘的总使用电力。

构成为：由该电力计5测量出的与总使用电力有关的信号被发送到工作估计部6，此处进行哪个电设备正在工作之类的关于电设备的机型的确定，并对其工作状态进行推测，将其推测结果显示在显示器7上。电力计5、工作估计部6以及显示器7构成本实施例的个体电设备工作估计装置。

上述工作估计部6例如由微型计算机构成，在图2中示出了其结构。

即，如图2所示，工作估计部6具备时间变动计算部6a、电设备机型确定部6b、对应似然度估计部6c、个体电设备工作状态估计部6d以及电设备的列表6g。

下面，分别对它们进行说明。

时间变动计算部6a被输入由电力计5测量出的与房屋1中的总使用电力有关的信号，针对该信号以规定时间(在本实施例中为1分钟)为间隔计算当前的总使用电力与前一个总使用电力之间的时间变动量、即总使用电力的时间变动1级差，并将其输出到工作电设备机型确定部6b、对应似然度估计部6c以及个体电设备工作状态估计部6d。

此外，时间变动计算部6a相当于本发明的时间变动计算单元。

电设备机型确定部6b基于由时间变动计算部6a计算出的时间变动1级差(相当于跳跃)来估计个体电设备的标准电力，对发生频度高的变动量(跳跃电力)进行聚类来形成簇，估计每个簇(机型)的跳跃电力的均值来作为各个个体电设备的标准电力。

在此，在新机型的判定基准中，在跳跃量超过与簇对应的标准电力的方差的某个固定比率的情况下判定为新机型，否则判定为现有的机型的电力使用处于标准电力的通常的变动范围内。

即，电设备机型确定部6b提取是否存在新的电设备，在提取出新的电设备的情况下作为新的电设备的机型追加到现有的机型列表6g，在未提取出新机型时维持现有的机型列表。

另外，在提取出新的电设备的情况下，按个体电设备的机型重新计算标准电力，根据需要来从机型列表6g删除(删减)实质上重复的机型等多余的机型。在此，机型的簇数是可变的。实时地逐渐学习与跳跃相应的个体电设备的标准电力的估计。将这样得到的新的机型列表6g的数据分别输出到对应似然度估计部6c、个体电设备工作状态估计部6d以及机型列表6g。此外，电设备机型确定部6b相当于本发明的电设备机型确定单元。

对应似然度估计部6c针对由电设备机型确定部6b得到的最新的机型列表的每个机型，估计关于由时间变动计算部6a计算出的时间变动1级差、即升(测量出的跳跃的符号为正)和降(测量出的跳跃的符号为负)这样的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度(概率)。该似然度被输出到个体电设备工作状态估计部6d。

此外，对应似然度估计部6c相当于本发明的对应似然度估计单元。

个体电设备工作状态估计部6d基于由时间变动计算部6a计算出的时间变动1级差(跳跃)以及由对应似然度估计部6c估计出的似然度，来实时地逐渐学习哪个电设备正在工作、以及其工作状况。为此，个体电设备工作状态估计部6d具备马尔可夫转换模型(Markovswitching model)6e和卡尔曼滤波器(Kalman filter)6f。

即，使用前者，在输入了在当前的个体电设备的工作状态下得到总使用电力的变动跳跃的、与个体电设备的启动/关闭的事件发生有关的似然度(概率)时估计什么样的个体电设备的工作状态发生变化，使用后者来根据总使用电力的时间变动的测定数据和上述工作状态的变化来优化个体电设备的机型估计、标准电力估计以及工作状态的估计。然后，按工作日/周六日、按时段对个体电设备的工作状态的估计结果进行统计，来学习个体电设备的使用模式。为了根据个体电设备的工作状态的估计结果及其过去的模式来逐渐修正先验分布，将该学习结果输出到对应似然度估计部6c、显示器7。

此外，个体电设备工作状态估计部6d相当于本发明的个体电设备工作状态估计单元。

电设备的列表6g存储所提取的机型的不同的标准电力的值(均值、方差)。在此存储的数据被发送到电设备机型确定部6b来用于此处的运算，接收运算的结果并根据需要重写数据。

接着，基于图3的流程图来说明由上述工作状态估计部6执行的个体电设备的工作状态估计处理的流程。

在个体电设备的工作状态估计中，如上所述那样使用NIALM实时逆估计算法。在此，作为表现出总使用电力的时间变动的动态估计模型的核心的部分一般是被称为状态空间模型(State Space Model)的模型。

该逆估计算法的最终输出是：以与测定出的总使用电力的时间变动的输入流对应的时间间隔(在本实施例中为1分钟)，来估计时刻变化的个体电设备的启动和关闭的工作状态。

像这样着眼于作为总使用电力的时间变动的跳跃是基于以下考虑：在电设备被启动时发生标准电力的垂直跳跃(符号为正的跳跃)，当在某一个时间点电设备被关闭时，发生抵消垂直跳跃的、符号为负的负跳跃，电设备的使用电力变为0。此外，这已如图9所示那样经实际的事例所验证。

在此，需要事先明确地定义作为状态变量的个体电设备的工作状态。在此，进行属于所估计出的个体电设备的“机型”的电设备的台数为1台这样的简单化的假设，将个体电设备的工作状态的状态变量定义为其工作概率，来进行模型化。然而，与各机型对应的个体电设备的台数为1台这样的简单化的假设是能够容易地扩展的，这是不言而喻的。

下面，按照图3所示的流程图来进行说明，首先在步骤S1中，工作状态估计部6从电力计5接收房屋1的总使用电力的数据。接着，进入步骤S2。

在步骤S2中，时间变动计算部6a基于通过步骤S1得到的总使用电力来计算以规定时间(在此为1分钟)为间隔的总使用电力的时间变动1级差的值(跳跃电力)，将该跳跃电力的数据分别输出到工作电设备机型确定部6b、对应似然度估计部6c以及工作状态估计部6d。接着，进入步骤S3。

在步骤S3中，工作电设备机型确定部6b使用从最初起逐渐学习并存储的机型(相当于簇)的列表6g、各机型的标准电力(其均值和方差)、工作概率的数据，基于通过步骤S2得到的跳跃电力来判定是否为与到此为止存储的机型不同的新机型。在该判定结果判定为是到此为止没有的新机型(判定为“是”)的情况下，进入步骤S4，在该判定结果判定为不是新机型(判定为“否”)的情况下，进入步骤S5。

具体地说，在最初的时间点机型数少至0种至1种的情况下，放宽追加新机型的判定基准，在经过时间长而能够取得足够的数据长度的时间点，使用信息量基准等严格的统计性判定基准。在经过时间短的时间点，例如，若是标准电力的大小的10％以下的跳跃电力则仍设为同一机型列表6g，而针对大于标准电力的大小的10％的跳跃电力则作为新机型追加到机型列表6g。

在步骤S4中，增加所存储的现有的机型(簇)的数量k。将通过步骤S2得到的跳跃电力的绝对值决定为新机型的标准电力(相当于簇的均值)，重新更新机型列表6g。接着，进入步骤S5。

在步骤S5中，与新增加的簇数k相应地重新计算每个机型的标准电力的均值。即，将总使用电力的以1分钟为间隔的时间序列数据中的、已决定的过去的数据长度之前的数据删除。具体地说，例如删除过去1天的1,440分钟或者1周的10,080分钟之前的数据，并且提取跳跃电力不是0的数据来作为要利用的数据。使用所提取出的该数据来重新计算向机型列表6g追加的机型的标准电力(簇的均值)。接着，进入步骤S6。

在步骤S6中，基于通过步骤S5重新计算的标准电力的均值的计算结果，来进行簇的重新分析，根据是否为同一标准电力、即是否为相同机型来判定是否对重复的多余机型进行删减。如果该判定结果是“是”，则进入步骤S7，如果判定结果是“否”，则进入步骤S8。

在步骤S7中，在从现有的机型列表删除了重复的多余机型的情况下，按删除了多余机型后的机型列表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测定出的总使用电力进行比较，实时地优化个体电设备的机型估计、标准电力估计以及工作状态的估计。

如以上那样，个体电设备的工作状态估计处理大致分为以图3所示的区域A、B、C示出的3个处理。

即，该图中的区域A包括步骤S2～S9，对应于如下处理：根据总使用电力的时间变动1级差数据来提取并确定工作电设备的机型。该区域A的处理是为了以下目的而执行的：判定是否将新的个体电设备新追加到现有的电设备的机型列表，另外，从现有的机型列表删除被认为是相同机型的机型。

另一方面，区域B包括步骤S9，对应于如下处理：对关于总使用电力的时间变动1级差、即升和降的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度进行估计。

另外，区域C包括各步骤S9～S11，对应于如下处理：以测定出的总使用电力的时间变动的时间间隔来动态地确定时刻变化的个体电设备的工作状态。

进一步详细地说，区域A的处理包括以区域A-1示出的步骤S3～S4的用于提取新机型的处理、以及步骤S4的标准电力的重新计算和多余机型的删减的处理。

另外，区域C的处理包括步骤S10的在输入了在当前的个体电设备的工作状态下根据总使用电力的跳跃而得到的、与个体电设备的启动/关闭的事件发生有关的似然度(概率)时使用马尔可夫转换模型6e来估计个体电设备的工作状态如何变化的处理、以及步骤S11的使用卡尔曼滤波器6f来根据总使用电力的时间变动的测定数据优化个体电设备的机型估计、标准电力估计及工作状态的估计的处理。

在图3所示的上述个体电设备的工作状态估计处理中，该整体的流程遵循状态空间模型的建模方法。

即，以上述区域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算法。

此外，关于下面的符号，为了方便而使说明文章中的符号不同于[数1]～[数113]中的符号和字体，但是它们是相同的。另外，在符号的上标和下标上下重叠的情况下，在说明文中左右(使下标在前，之后为上标)错开地记载。

从混合正态模型视点而言，在本实施例的情况下k-均值模型是单纯的1维模型。

即，将m个样本表示为x_i,(i＝1，…m)，另外，将k个未知的正态分布表示为N(μ_j，σ_j ²)，(j＝1，…，k)。在此k已知。此外，为了便于说明，在此设σ_j,(j＝1,…,k)是已知的，σ_j＝σ^{^},(j＝1，…，k)且全部相等。在图4中示出了k＝2的情况下的示意图。

根据该图可以明确，从样本来看，x_i从2个正态分布N(μ_j,σ^{^2}),j＝1,2得到的似然度分别如下。

[数1]

$l_{ij}=\frac{1}{\sqrt{2\mathrm{\π}}\widehat{\mathrm{\σ}}}{e}^{-\frac{1}{2{\widehat{\mathrm{\σ}}}^2}{(x_i-{\mathrm{\μ}}_j)}^2},j=1,2...\left(1\right),$

现在，若将收敛计算的本期第步的均值μ_j,j＝1,2的估计值设为μ_j ^(s),j＝1,2，则可将式(1)认作在第s步中从2个正态分布N(μ_j,σ^{^2}),j＝1,2得到样本x_i的似然度的估计值，因此将其同样地表示为l_y ^(x),(j＝1,2；i＝1,…,m)。在利用k-均值模型的EM算法的解法中，通过下面的加权平均的收敛计算来求出使各样本从未知的2个正态分布得到的似然度最大化的未知参数μ_j,j＝1,2。

[数2]

${\mathrm{\μ}}_j^{(s+1)}=\frac{\underset{i=1}{\overset{m}{\Σ}}{l}_{ij}^{\left(s\right)}{x}_i}{\underset{i=1}{\overset{m}{\Σ}}{l}_{ij}^{\left(s\right)}},j=1,2;s=0,1,......\left(2\right)$

根据式(1)，式(2)的右边是μ_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的值的潜在随机变量。

[数3]

那么，设想一个不可观测的随机变量Z_ij，下面给出基于假设观测到随机变量Z_ij时的完整数据y_i,i＝1,…,m的表达。

完整数据：

[数4]

y_i＝(x_i，z_ij，j＝l，...，k)，i＝1，...，m ...(4)

统计模型：

[数5]

$x_i~N({\mathrm{\μ}}_j,{\mathrm{\σ}}_j^2),j=1,...,k...\left(5\right)$

[数6]

如上所述，在此设样本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]时，该假设下的完整数据的似然度为下式。

[数7]

$\begin{array}{c}p\left(y_i\right|h^{\′})=\underset{j=1}{\overset{k}{\mathrm{\Π}}}{\left(\frac{1}{\sqrt{2{\mathrm{\π\σ}}_j^{\′2}}}{e}^{-\frac{1}{2{\mathrm{\σ}}_j^{\′2}}{(x-{\mathrm{\μ}}_j^{\′})}^2}\right)}^{Z_{ij}}\\ =\underset{j=1}{\overset{k}{\mathrm{\Π}}}{\left(\frac{1}{\sqrt{2{\mathrm{\π\σ}}_j^{\′2}}}\right)}^{Z_{ij}}{e}^{-\frac{1}{2}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{Z}_{ij}\frac{{(x-{\mathrm{\μ}}_j^{\′})}^2}{{\mathrm{\σ}}_j^{\′2}}}\end{array}...\left(7\right)$

若为了简化计算而取两边的对数来求出全部样本的似然度，则为下式。其中，在此y＝(y_i,i＝1,…,m)为完整数据矢量。

[数8]

$\begin{array}{c}logp\left(y\right|h^{\′})=log\underset{i=1}{\overset{m}{\Π}}p\left(y_i\right|h^{\′})\\ =\underset{i=1}{\overset{m}{\Π}}logp\left(y_i\right|h^{\′})\\ =\underset{j=1}{\overset{k}{\mathrm{\Σ}}}(-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{z}_{ij}\log (2{\mathrm{\π\σ}}_j^{\′2}))-\frac{1}{2}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\Σ}}{z}_{ij}\frac{{(x-{\mathrm{\μ}}_j^{\′})}^2}{{\mathrm{\σ}}_j^{\′2}}\\ =\underset{j=1}{\overset{k}{\mathrm{\Σ}}}(-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{z}_{ij}\log 2\mathrm{\π}-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{z}_{ij}{\mathrm{log\σ}}_j^{\′2}-\underset{i=1}{\overset{m}{\Σ}}\frac{1}{2}{z}_{ij}\frac{{(x_i-{\mathrm{\μ}}_k^{\′})}^2}{{\mathrm{\σ}}_j^{\′2}})\end{array}...\left(8\right)$

使该对数似然度最大化的、(μ_j,σ_j ²),j＝1,…,k的最大似然估计值(μ_j ^{^}’,σ_j ^{^，2}),j＝1,…,k是通过使式(8)关于(μ_j’,σ_j’²),j＝1,…,k最大化而求出的。分别如下。

关于均值μ_j’,j＝1,…,k的最大似然估计值μ_j ^{^}’,j＝1,…,k，如果对式(8)针对μ_j’,j＝1,…,k进行微分并根据一阶条件求出μ_j ^{^}’,j＝1,…,k，则得到下式。

[数9]

${{\widehat{\mathrm{\μ}}}^{\′}}_j=\frac{\underset{i=1}{\overset{m}{\Σ}}{Z}_{ij}{x}_j}{\underset{i=1}{\overset{m}{\Σ}}{Z}_{ij}},j=1,...,k...\left(9\right)$

另外，关于方差σ_j’²,j＝1,…,k的最大似然估计值σ_j ^{^，2},j＝1,…,k，同样地，当以σ_j’²,j＝1,…,k对式(8)进行微分并求出一阶条件时，得到下式。

[数10]

$\frac{\underset{i}{\Σ}{Z}_{ij}}{2}\·\frac{1}{{\mathrm{\σ}}_j^{\′2}}\·(-1)\frac{1}{2}\·(-1)\frac{1}{{\left({\mathrm{\σ}}_j^{\′2}\right)}^2}\underset{i=1}{\overset{m}{\Σ}}{Z}_{ij}{(x_i-{\mathrm{\μ}}_j^{\′})}^2=0...\left(10\right)$

若对其进行求解则得到下式。

[数11]

${\widehat{\mathrm{\σ}}}_j^{\′2}=\frac{\underset{i=1}{\overset{m}{\Σ}}{Z}_{ij}{(x_i-{{\widehat{\mathrm{\μ}}}^{\′}}_j)}^2}{\underset{i=1}{\Σ}{z}_{ij}}...\left(11\right)$

因而，如果将最大似然估计值(μ_j ^{^}’,σ_j ^{^，2}),j＝1,…,k代入式(8)则能够求出最大对数似然度l*。

[数12]

在此，当代入

[数13]

${\widehat{\mathrm{\σ}}}_j^{\′2}=1/\left({\mathrm{\Σ}}_{i=1}^m{Z}_{ij}\right)\·{\mathrm{\Σ}}_{i=1}^m{Z}_{ij}{(x_i-{{\widehat{\mathrm{\μ}}}^{\′}}_j)}^2$

时，能够如下那样表示式(12)。

[数14]

那么，在本实施例中，考虑了簇的数量k不固定而可变的情况。在决定最合适的簇数k*时，能够利用赤池的信息量准则(AIC)。

能够使用该AIC来将最合适的簇数设为使AIC最小的k*。当将最大对数似然度表示为l*、将参数数量表示为p时，能够用下式来定义AIC。

[数15]

需要注意的是，在k-均值模型中，当使簇的数量k增加一个时会增加均值和方差这2个参数。另外，最大似然估计值(μ_j ^{^}’,σ_j ^{^，2}),j＝1,…,k、最大对数似然度i*、以及AIC也都依赖于簇的数量k而改变，因此将它们如AIC(k)那样表示为k的函数。

作为求出最合适的簇数k*的算法的表达，只要依次使k从1逐渐增加，将满足AIC(k+1)>AIC(k)的k设为最合适的簇数k*即可。

具体地说，当求出AIC(k+1)、AIC(k)时，得到下式(15)、(16)。

[数16]

[数17]

因而，相当于求出使下式成立的k。

[数18]

在不完整数据最大似然估计法中，将EM算法或将其扩展为点对集映射而得到的GEM(Generalized Expectation Maximization：广义期望最大化)算法应用于不完整数据下的重力模型的估计、IPFP(Iterative Proportional Fitting Procedure：迭代比例拟合程序)的扩展。首先，说明不完整数据的定义、GEM算法的特征性。

考虑随机变量Y和X。数据是指这些随机变量的实现值，以小字符y、x表示。一般来说，在随机变量y与x之间存在基于多对一的映射h的对应关系时，将y定义为完整数据，将x定义为不完整数据。即，

[数19]

X＝h(Y) 多对一 …(18)

当利用数据来表达时同样，在完整数据y与不完整数据x之间设想以下的关系。

[数20]

x＝h(y) 多对一 …(19)

该定义是一般性的抽象的定义，但是为了可靠地理解，事先确认：潜在随机变量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的组，

[数21]

y＝(x，z) …(20)

。在此，如下式那样定义多对一的映射h。

[数22]

x＝h(x，z)≡h(y) …(21)

实际上，无论针对z所取的任何值，只要确定了与z成对的x的值，映射h的像就始终为x，据此可知映射h为多对一。关于这一点，若想像将(x,z)的2维平面投影到x轴上的情况则易于理解。

那么，在可观测的只有随机变量x时、换言之在可观测的只有不完整数据x时，在最大似然估计法中需要求出使生成该不完整数据的似然度最大的未知参数的值来作为该参数的最大似然估计值。

该方法如下：在某个模型(参数)下，将生成不完整数据x的概率(似然度)返回到完整数据y的模型来进行计算。

也就是说，若注意到与映射的像x对应的原像是生成像x的完整数据y的集合，则生成不完整数据x的概率(似然度)是原像h^-1(x)所包含的完整数据y发生的概率。换言之，原像h^-1(x)所包含的各个完整数据y发生的概率(似然度)之和为不完整数据x发生的似然度。当通过数式来表达以上内容时如下。

首先，下面定义使不完整数据x发生的完整数据y的集合。

[数23]

h^-1(x)＝(y|x＝h(y)} …(22)

若将不完整数据x发生的概率(似然度)表示为g(x|Φ)则如下。

[数24]

$g\left(x\right|\mathrm{\φ})={\∫}_{h^{-1}\left(x\right)}f\left(y\right|\mathrm{\φ})dy...\left(23\right)$

因而，不完整数据的最大似然估计法是求出使式(23)的似然度最大的参数Φ。

一般来说，使不完整数据的似然度g(x|Φ)关于Φ最大化的问题的解法并不简单。作为该解法，考虑了被称为EM(Expectation Maximization：期望最大化)、GEM(Generalized Expectation Maximization：广义期望最大化)的方法。

这些方法的概要如下。

首先，以L(Φ)来表示不完整数据的对数似然度。

[数25]

L(φ)＝log(g(x|φ)) …(24)

接着，以k(y|x,Φ)来表示当给定不完整数据x和参数Φ时的完整数据y的条件概率。

[数26]

k(y|x，φ)＝f(y|φ)/g(x|φ) …(25)

另外，分别如下那样定义Q(Φ’|Φ)和H(Φ’|Φ)。

[数27]

Q(φ′|φ)＝E[log(f(y|φ′))|x，φ] …(26)

[数28]

H(φ′|φ)＝E[log(k(y|x，φ′))|x，φ] …(27)

根据以上的定义，用下式表示对数似然度L(Φ’)。

[数29]

L(φ′)＝Q(φ′|φ)-H(φ′|φ) …(28)

这是由于，一般来说，E[h(x)|x]＝h(x)成立，由此

[数30]

$\begin{array}{c}L\left({\mathrm{\φ}}^{\′}\right)=log\left(g\right(x\left|{\mathrm{\φ}}^{\′}\right))=E\[log\left(g\right(x\left|{\mathrm{\φ}}^{\′}\right))|x,\mathrm{\φ}\]\\ =E\[log\left(f\right(y\left|{\mathrm{\φ}}^{\′}\right)\left)\right|x,\mathrm{\φ}\]-E\[log\left(k\right(y|x,{\mathrm{\φ}}^{\′})\left)\right|x,\mathrm{\φ}\]\end{array}...\left(29\right)$

成立。

EM算法是包括下面两个步骤的、根据第步中的Φ的估计值Φ⁽ⁿ⁾来计算第步的估计值Φ⁽ⁿ⁺¹⁾的Φ⁽ⁿ⁾→Φ⁽ⁿ⁺¹⁾,n＝0,1,2,…的迭代计算法。在此，在判定是否为新的步的期望值步骤(E-步骤)和对均值进行评价的最大化步骤(M-步骤)中，分别进行以下的计算。

EM算法的定义：

E-步骤：

[数31]

Q(φ|φ⁽ⁿ⁾) …(30)

M-步骤：

[数32]

${\mathrm{\φ}}^{(n+1)}=\underset{\mathrm{\φ}\∈\mathrm{\Ω}}{\max }Q\left(\mathrm{\φ}\right|{\mathrm{\φ}}^{\left(n\right)})...\left(31\right)$

另一方面，GEM算法是将EM算法的M步骤置换为下面的点对集映射M而得到的。

GEM的M-步骤：

[数33]

GEM的M-步骤：φ⁽ⁿ⁾→φ⁽ⁿ⁺¹⁾∈M(φ⁽ⁿ⁾) …(32)

在此，如下那样定义点对集映射M。

[数34]

M(φ)＝{φ′∈Ω|Q(φ′|φ)≥Q(φ|φ)} …(33)

此外，已知下式成立，

[数35]

H(φ′|φ)≤H(φ|φ) …(34)

使用该式来证明了GEM的单调(增加)性。

[数36]

L(M(φ))≥L(φ) …(35)

因而，如果对数似然度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(Φ’|Φ)。能够参照式(8)来如下那样表示Q(h’|h)。

[数37]

$\begin{array}{c}Q\left(h^{\′}\right|h)=E\[logp\left(y\right|h^{\′})|h\]\\ =\underset{z}{E}\[logp(x,z|z,{\mathrm{\μ}}^{\′},{\mathrm{\σ}}^{\′2})|x,\mathrm{\μ},{\mathrm{\σ}}^2\]\\ =E\[-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{z}_{ij}\log 2\mathrm{\π}-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{z}_{ij}{{\mathrm{log\σ}}^{\′}}_j^2-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{z}_{ij}\frac{{(x_i-{\mathrm{\μ}}_j^{\′})}^2}{{\mathrm{\σ}}_j^{\′2}}|x,\mathrm{\μ},{\mathrm{\σ}}^2\]\\ =-\frac{1}{2}m\log 2\mathrm{\π}-\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}E\left(z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2){{\mathrm{log\σ}}^{\′}}_j^2\\ -\frac{1}{2}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\Σ}}E\left(z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)\frac{{(x_i-{\mathrm{\μ}}_j^{\′})}^2}{{\mathrm{\σ}}_j^{\′2}}\end{array}...\left(36\right)$

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步骤中，下面的计算成为问题。

[数38]

$E\left(Z_{ij}\right|X_i=x_i,\mathrm{\μ}=\left({\mathrm{\μ}}_1,...,{\mathrm{\μ}}_k\right),{\mathrm{\σ}}^2=\left({\mathrm{\σ}}_1^2,...,{\mathrm{\σ}}_k^2\right))...\left(37\right)$

[数39]

在此，说明式(37)的计算方式。

式(37)采用后验概率的形式。即，在观测到不完整数据x_i的基础上，求出指示器随机变量Z_g的后验概率，该指示器随机变量Z_g用于判别x_i是来自哪个正态分布的样本。虽然不知道该后验概率，但是在给定指示器随机变量Z_g的基础上，能够明确地计算不完整数据x_i发生的似然度。

也就是说，通过使用贝叶斯定理，设式(37)的后验概率同先验分布与似然度之积成正比的值来求出式(37)的后验概率。

因此，当将Z_g＝1的事件表示为Ω_y时，下式成立。

[数40]

Ω_ij＝{ω∈Ω|Z_ij＝1} …(39)

另一方面，对于全部的j≠j’，Ω_y与Ω_y’明显是互斥的，因此，

[数41]

另外，

[数42]

$\underset{j=1}{\overset{k}{\∪}}{\mathrm{\Ω}}_{ij}=\mathrm{\Ω}...\left(41\right)$

成立。因而，在期望值与概率之间成立下面的关系。

[数43]

E(Z_ij|x，μ，σ²)＝1·P(Z_ij＝1|X＝x_i，μ，σ²) ···(42)

当使用该式时，根据贝叶斯定理，得到下式。

[数44]

P(Z_ij|X_i)∝P(X_i|Z_ij)P(Z_ij) …(43)

假定按照时间顺序t观测到不完整数据x_i，来将先验概率P(Z_y)扩展为依赖于紧挨着t期之前的X_t-1,…,X₁的观测的先验分布。

当设为πj|i-1＝P(Z_y|X_i-1,…,X₁)时，式(43)表示为下式。

[数45]

$E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)=\frac{\frac{1}{\sqrt{2{\mathrm{\π\σ}}_j^2}}{e}^{-\frac{1}{2}\·\frac{{(x_i-{\mathrm{\μ}}_j)}^2}{{\mathrm{\σ}}_j^2}}\·{\mathrm{\π}}_{j|i-1}}{\underset{j=1}{\overset{k}{\Σ}}\frac{1}{\sqrt{2{\mathrm{\π\σ}}_j^2}}{e}^{-\frac{1}{2}\·\frac{{(x_i-{\mathrm{\μ}}_j)}^2}{{\mathrm{\σ}}_j^2}}\·{\mathrm{\π}}_{j|i-1}}...\left(44\right)$

在设为

[数46]

${\mathrm{\φ}}_j\left(x_i\right)=1/\left(\sqrt{2{\mathrm{\π\σ}}_j^2}\right)\exp \left(\right(-1/2)\·({\left(x_i-{\mathrm{\μ}}_j\right)}^2/{\mathrm{\σ}}_j^2\left)\right)$

的情况下，能够如下那样总结作为本实施例中使用的k-均值的解法算法的EM算法。

E步骤：

[数47]

$E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)=\frac{{\mathrm{\φ}}_j\left(x_i\right){\mathrm{\π}}_{j|i-1}}{\underset{j=1}{\overset{k}{\Σ}}{\mathrm{\φ}}_j\left(x_i\right){\mathrm{\π}}_{j|i-1}}...\left(45\right)$

M步骤：

[数48]

${\widehat{\mathrm{\μ}}}_j^{\′}=\frac{\underset{i=1}{\overset{m}{\Σ}}E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)x_i}{\underset{i=1}{\overset{m}{\Σ}}E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)}...\left(46\right)$

[数49]

${\widehat{\mathrm{\σ}}}_j^{\′2}=\frac{\underset{i=1}{\overset{m}{\Σ}}E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2){(x_i-{\widehat{\mathrm{\μ}}}_j^{\′})}^2}{\underset{j=1}{\overset{m}{\Σ}}E\left(Z_{ij}\right|x,\mathrm{\μ},{\mathrm{\σ}}^2)}...\left(47\right)$

根据式(46)和式(47)可知，该EM算法为以h＝[μ,σ²]为输入来生成h’＝[μ^{^}’,σ^{^，2}]的迭代计算法。

在以上的说明中，为了避免复杂，考虑了固定了簇的数量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的基准进行个体电设备的机型的增加、删除等，则会立刻失败。该问题一般也被称为冷启动(coldstart)的问题，作为以下问题而为人所知：即使是推荐系统(recommendation system)等的系统设计，对于新用户来说，利用历史记录少，从而也发生系统工作不顺利的情况。

因此，在本实施例中，考虑到冷启动问题，在簇数的决定中采纳用于解决该问题的实用的方法。

在此，簇数表示个体电设备的机型的数量，另外，构成簇的正态分布的均值与个体电设备的标准电力对应。下面，示出本实施例中采纳的可变簇数的决定方法。

在本实施例中，首先着眼于总使用电力的以1分钟为时间间隔的级差的电力跳跃量。

图3的区域A-1中的、根据跳跃量来识别新的电设备的机型并追加到现有的机型列表的处理中，

1)在跳跃量接近0的情况下，

认为个体电设备的工作状态未发生任何变化，转移到下一期。

另一方面，

2)在跳跃不是0的情况下，

a)在跳跃量接近现有的电设备的标准电力的情况下，判定为现有的电设备被启动或关闭，

b)在跳跃偏离于现有的电设备的标准电力的情况下，判定为确定出新的电设备的机型，将新的电设备的机型的标准电力设定为该跳跃量并将伴有该跳跃量的新机型追加到现有的机型列表。

这样，在上述的新的电机型的导入判定处理中，是1)对现有的电设备的机型列表追加新机型、或者2)什么都不追加中的任一个处理。即，作为判定为新机型的基准，在跳跃量超过了标准电力的方差的某个固定比率的情况下判定为新机型，否则判定为现有的机型的电力使用处于标准电力的通常的变动范围内，从而不变更现有的机型列表。

基于通过区域A-1的处理而决定的电设备的机型列表来执行EM算法，进行个体电设备的标准电力的更新、重新计算。

这是下面的删除(删减)电机型的重复的处理。

即，在区域A-2的处理中，更新标准电力及其方差以及删除机型的重复。

1)基于通过区域A-1的处理而得到的电设备的机型列表及其标准电力以及标准电力的方差使EM算法进行1步计算，由此进行个体电设备的标准电力和标准电力的方差的基于重新计算的更新。

2)根据标准电力的重新计算的结果，将视为标准电力相同的机型合并为一个机型。即，删除重复机型。

这样，区域A-2的处理是进行利用EM算法的均值、方差的重新计算、更新的处理。在此也同样地，关于是否属于同一机型的判定，以与标准电力的方差相对的大小为判定基准。

那么，EM算法原本是进行几次迭代收敛计算并将收敛得到的值作为均值、方差的估计值的算法。但是，如果使用以1分钟为间隔而输入的总使用电力、其跳跃量来每隔1分钟地反复进行收敛计算直到各期中收敛为止，则计算成本非常高。

因此，在本实施例中，利用具有似然度的单调性的特征的EM算法的特性，将EM算法的迭代计算仅进行1次，并继承到下一个处理。

如果接近同一个此前登记的标准电力的值，则认为与该标准电力对应的电设备被启动、关闭。是否接近现有的标准电力的值的判定是利用与标准电力对应的方差的大小来进行判定的。

此外，在上述M步骤中，将所使用的过去的数据的长度、返回的长度指定为比所述过去的数据少的长度。

以上，结束了簇数可变的可变k-均值算法的说明，因此接着说明用于NIALM的逆估计算法的动态表达。

用于NIALM的实时逆估计算法的目的在于，在实际的家庭环境中，根据总使用电力的以1分钟为时间间隔的变动数据(跳跃量及其升方向或降方向)，实时地估计存在什么样的电设备、处于什么样的工作状态。为此，需要实时且动态地处理以时间序列变动的数据的系统。并且，不仅需要实时地处理时间序列数据，还需要结合实时的控制。

目前，作为一边实时且动态地处理时间序列数据、一边结合控制的方法，作为最有力的方法，已知一种被称为状态空间建模的方法。

使用该状态空间模型的方法也被称为卡尔曼滤波器。其基本考虑方法包括：表达作为潜在变量的状态变量的时间性转移的状态方程；以及将不可观测的状态变量与观测数据相结合的观测方程。

在本实施例中，状态方程是与根据总使用电力的1级级差的跳跃来估计哪个电设备被启动/关闭之类的统计性估计有关的模型的描述式。

考虑状态空间建模时的重要事项是：1)如何定义状态变量；2)描述状态如何随时间而转移的状态方程；以及3)将作为不可观测的潜在变量的状态变量与观测变量相结合的观测方程。

并且，4)在与控制相关联的情况下，除了状态变量以外，控制者、观测者需要对以状态方程表示的动态系统追加输入政策变量或控制变量。

本实施例的用于NIALM的实时逆估计算法的特征在于进行包含该最后的输入变量的模型化。即，在本实施例中的包含输入变量的模型化中，使用不是以往的控制变量的意思、而是作为改变动态系统的状态的外在冲击的模型化。若更广义地一般化，则是对来自外部的学习输入等也能够进行同样的处理的方法。

在此，示出卡尔曼滤波器的表达的标准形。

首先，将t期的状态变量矢量设为x_i，将输入矢量设为u_t，将观测矢量设为y_t。将状态方程的误差项的随机变量矢量设为W_t，将观测方程的误差项的随机变量矢量表示为v_t。

在标准形中，上述的状态方程、观测方程分别如下那样表示。

[数50]

x_t＝F_t-1x_t-1+G_t-1u_t-1+w_t-1 …(48)

[数51]

y_t＝H_tx_t+v_t …(49)

其中，F_t、G_t、H_t是已知的系数矩阵。

在此，状态方程的误差项以及观测方程的误差项在不同时间点之间是独立的。另外，状态方程的误差项与观测方程的误差项彼此独立。

下面表达以上的假定。

[数52]

$E(w_t,w_{t^{\′}}^T)=Q_t{\mathrm{\δ}}_{{\mathrm{tt}}^{\′}}...\left(50\right)$

[数53]

$E(v_t,v_{t^{\′}}^T)=R_t{\mathrm{\δ}}_{{\mathrm{tt}}^{\′}}...\left(51\right)$

[数54]

$E(w_t,v_{t^{\′}}^T)=O...\left(52\right)$

其中，δ_n’是在t＝t’时为1、除此以外时为0的克罗内克函数(Kronecker delta)。另外，式(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)没有区别。

在此，引入下面的条件概率的表示法。这样一来，能够使卡尔曼滤波器的递归修正的表达简明化。

[数55]

${\hat{x}}_t^{-}=E\left(x_t\right|y_1,y_2,...,y_{t-1})...\left(53\right)$

[数56]

${\hat{x}}_t^{+}=E\left(x_t\right|y_1,y_2,...,y_{t-1},y_t)...\left(54\right)$

即，在上述各式中，右上标表示在给定1期前为止的观测值时的条件期望值，右上标表示在给定本期为止的观测值时的条件期望值。

同样地，状态变量矢量的协方差矩阵的条件期望值中也使用同样的表示法。

[数57]

$P_t^{-}=E\&lsqb;(x-{\hat{x}}_t^{-}){(x_t-{\hat{x}}_t^{-})}^T|y_1,y_2,...,y_{t-1}\&rsqb;...\left(55\right)$

[数58]

$P_t^{+}=E\&lsqb;(x_t-{\hat{x}}_t^{+}){(x_t-{\hat{x}}_t^{+})}^T|y_1,y_2,...,y_t\&rsqb;...\left(56\right)$

将其简记来如下那样表示。

[数59]

$P_t^{-}=E\&lsqb;(x_t-{\hat{x}}_t^{-}){(x_t-{\hat{x}}_t^{-})}^T\&rsqb;...\left(57\right)$

[数60]

$P_t^{+}=E\&lsqb;(x_t-{\hat{x}}_t^{+}){(x_t-{\hat{x}}_t^{+})}^T\&rsqb;...\left(58\right)$

P_t ^-、P_t-⁺均为状态变量矢量的误差协方差。

如下那样给出初始步骤。

[数61]

${\hat{x}}_0^{+}=E\left(x_0\right)...\left(59\right)$

[数62]

$P_0^{+}=E\&lsqb;(x_0-{\hat{x}}_0^{+}){(x_0-{\hat{x}}_0^{+})}^T\&rsqb;...\left(60\right)$

在此，输入的序列u₀,u₁,…u_t,…是问题，假定输入的序列是无误差地从外部给定的。

作为卡尔曼滤波器的核心的递归性(Recursion)的更新处理如下那样表达。

[数63]

${\hat{x}}_t^{-}=F_{t-1}{\hat{x}}_{t-1}^{+}+G_{t-1}{u}_{t-1}...\left(61\right)$

[数64]

$P_t^{-}=F_{t-1}{P}_{t-1}^{+}{F}_{t-1}^T+Q_{t-1}...\left(62\right)$

[数65]

$K_t=P_t^{-}{H}_t^T{(H_t{P}_t^{-}{H}_t^T+R_t)}^{-1}...\left(63\right)$

[数66]

${\hat{x}}_t^{+}={\hat{x}}_t^{-}+K_t(y_t-H_t{\hat{x}}_t^{-})...\left(64\right)$

[数67]

$P_t^{+}=P_t^{-}-K_t{H}_t{P}_t^{-}...\left(65\right)$

式(63)是被称为卡尔曼增益的项。

遵循以上的方程体系的卡尔曼滤波器的计算过程如下。

1)首先，基于式(59)和式(60)来决定初始值(x₀ ^{^+},P₀ ⁺)。

2)接着，基于式(61)和式(62)，基于直到1期前为止的观测值来求出状态变量矢量的预测值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期的标准电力的估计值

(σ^{^2})_t ^k：第k种电设备的第t期的标准电力的方差的估计值

观测值是总使用电力的1维的以1分钟为间隔的观测数据。

y_t＝Pw_t,t＝1,2,…

此外，为了使记法与卡尔曼滤波器的标准形一致，使用y_t的符号。

那么，虽然是状态变量，但是由于用于NIALM的实时逆估计算法的目的是个体电设备的工作状态的动态估计，因此当然需要设状态变量矢量表示各电设备的工作状态。为了进一步明确这一点，如以下那样定义工作状态取启动/关闭的潜在随机变量。

将表示第种电设备的工作状态的随机变量设为X_t ^k。如以下那样定义状态随机变量X_t ^k。

[数68]

在此，当取状态随机变量的期望值X_t ^k时，成为

[数69]

$\begin{array}{c}{x}_t^k=E\left(X_t^k\right)=1\&CenterDot;P(X_t^k=1)+0\&CenterDot;P(X_t^k=0)\\ =P(X_t^k=1)\end{array}...\left(67\right)$

可知，期望值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的估计值的列表，即

[数70]

$B_t=\{(P_t^1,{\left({\mathrm{\&sigma;}}^2\right)}_t^1),(P_t^2,{\left({\mathrm{\&sigma;}}^2\right)}_t^2),...,(P_t^{k_t},{\left({\mathrm{\&sigma;}}^2\right)}_t^{k_t})\}$

因而，K_t＝|Ω_t|。其中，|S|表示集合S的元素的数量。

本实施例的逆估计算法的最大特征在于，将总使用电力的以1分钟为间隔的时间级差、即跳跃引入为卡尔曼滤波器的输入变量。

但是，并不是将跳跃直接作为输入，而是判定跳跃是基于现有的保有电设备的集合Ω_t中的哪个设备的工作状态的变化而产生的，使用k-均值模型来变换为针对各状态变量的输入矢量，这一点与卡尔曼滤波器的标准形不同。

首先，将作为总使用电力的1级级差的跳跃作为标量的输入变量，设为与总使用电力的状态变量X_t ⁰对应的标量的输入变量。

[数71]

$u_t^0=X_t^0-X_{t-1}^0...\left(68\right)$

通常，视为输入变量没有误差，因此u_t ⁰＝E(u_t ⁰)，下式成立。

[数72]

$u_t^0=E\left(X_t^0\right)-E\left(X_{t-1}^0\right)=x_t^0-x_{t-1}^0...\left(69\right)$

以后，对于输入变量，与期望值无区别地对待。

接着，定义与k_t种类的电设备对应的针对状态变量的输入变量u_t ^k,k＝1,…,k_t。

它们被定义为k_t个电设备的种类上的概率分布，该概率分布表示u_t ⁰的跳跃是基于1,…,k_t的哪个设备的工作状态的变化而产生。对于这些输入变量u_t ^k,k＝1,…,k_t，之后叙述形式上的定义，而关于其含义另行详细叙述。当前认为是与k_t个状态对应的输入矢量。

将这k_t+1个输入变量汇总起来表示为矢量u_t。

[数73]

$u_t=(u_t^0,u_t^1,...,u_t^{k_t})...\left(70\right)$

在以上的准备的基础上，构成状态方程。为此，关于x_t ^k,(k＝0,1,…,k_t)的k_t+1个状态变量，只要描述基于x_t-1 ^k和u_t ^k如何生成x_t ^k即可。

1)决定u_t ⁰和x_t ⁰

[数74]

$u_t^0={\mathrm{\&Delta;y}}_t=y_t-y_{t-1}={\mathrm{Pw}}_t-{\mathrm{Pw}}_{t-1}...\left(71\right)$

[数75]

$x_t^0=x_{t-1}^0+u_t^0...\left(72\right)$

根据该定义，状态变量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)。

[数76]

${\mathrm{nzu}}_t^k=\mathrm{\&phi;}\left(\right|u_t^0\left|\right|P_t^k,{\mathrm{\&sigma;}}_t^k){\mathrm{\&pi;}}_t^k/\underset{k=1}{\overset{k_t}{\mathrm{\&Sigma;}}}\mathrm{\&phi;}\left(\right|u_t^0|{|}_t{P}_t^k,{\mathrm{\&sigma;}}_t^k){\mathrm{\&pi;}}_t^k...\left(73\right)$

其中，

Φ(x|μ,σ)：均值μ、标准偏差σ的正态分布N(μ,σ²)的密度函数

|u_t ⁰|：u_t ⁰的绝对值

π_t ^k,k＝1,…,k_t：对k_t个假设N(μ_t ^k,(σ²)_t ^k)的先验分布(参见式(43))

式(73)表示判别总使用电力的跳跃是从哪个机型产生的概率。

(b)决定状态变量x_t ^k,k＝1,…,k_t(状态方程中的除去k＝0的、k＝1,…,k_t的方程)

下式示出了状态方程中除去k＝0的k＝1,…,k_t的方程。

[数77]

$\left\{\begin{array}{ccc}{u}_t^0>0,& x_t^k=x_{t-1}^k+\left(x_{t-1}^k\right){\mathrm{nzu}}_t^k+w_t^k& k=1,...,k_t\\ u_t^0<0,& x_t^k=x_{t-1}^k-x_{t-1}^k{\mathrm{nzu}}_t^k+w_t^k& k=1,...,k_t\\ u_t^0=0,& x_t^k=x_{t-1}^k+w_t^k& k=1,...,k_t\end{array}\right....\left(74\right)$

在此，用下式定义u_t ^k,(k＝1,…,k_t)。

[数78]

$\begin{array}{c}{u}_t^k=\left\{\begin{array}{ccc}(1-x_{t-1}^k){\mathrm{nzu}}_t^k& if& u_t^0>0\\ -x_{t-1}^k{\mathrm{nzu}}_t^k& if& u_t^0<0\\ 0& if& u_t^0=0\end{array}\right....\left(75\right)\\ k=1,...,k_t\end{array}$

下式回到卡尔曼滤波器的标准形。

[数79]

$\left\{\begin{array}{ccc}{u}_t^0>0,& x_t^k=x_{t-1}^k+u_t^k+w_t^k& k=1,...,k_t\\ u_t^0<0,& x_t^k=x_{t-1}^k+u_t^k+w_t^k& k=1,...,k_t\\ u_t^0=0,& x_t^k=x_{t-1}^k+w_t^k& k=1,...,k_t\end{array}\right....\left(76\right)$

3)决定观测方程

[数80]

$y_t=\underset{k=1}{\overset{k_t}{\mathrm{\&Sigma;}}}{P}_t^k{x}_t^k+v_t...\left(77\right)$

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的实时逆估计算法回到基于状态空间建模的卡尔曼滤波器的标准形。

更准确地说，若再次示出卡尔曼滤波器的标准形，则为下面的线性形式的标准形。

[数81]

x_t＝F_t-1x_t-1+G_t-1u_t+w_t …(78)

根据本项中的描述可以明确，并不是将总使用电力u_t ⁰直接代入到上面的式(48)，而是将经过式(73)的非线性的变换后的k_t个输入u_t ^k,(k＝1,…,k_t)代入到用式(48)来表达的k_t个状态转移方程，因此本质上是非线性的卡尔曼滤波器。

式(48)和式(76)是成为k-均值模型与卡尔曼滤波器的合并的核心的部分。式(73)是与图3中的区域BA的部分相当的部分，是计算如何将来自k-均值算法的输出变换为卡尔曼滤波器的输入的处理。

另一方面，式(74)是与图3的区域C-1的马尔可夫转换相当的部分，是计算根据上述输入如何改变状态变量的处理。

接着，详细说明这些计算处理。

首先，从将k-均值的输出变换为卡尔曼滤波器的输入的式(73)开始进行说明。

期中的k-均值模型的计算流程是以下流程：判定是否将电设备的新机型追加到现有列表，接着，基于该新的机型列表，通过EM算法来重新计算标准电力和标准电力的方差，如果有重复的机型则删除。

因而，来自图3的区域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。

[数82]

$\begin{array}{c}{\mathrm{nzu}}_t^k=\mathrm{\&phi;}\left(\right|u_t^0\left|\right|P_t^k,{\mathrm{\&sigma;}}_t^k){\mathrm{\&pi;}}_t^k/\underset{k=1}{\overset{k_t}{\mathrm{\&Sigma;}}}\mathrm{\&phi;}\left(\right|u_t^0\left|\right|P_t^k,{\mathrm{\&sigma;}}_t^k){\mathrm{\&pi;}}_t^k\\ k=1,...,k_t\end{array}...\left(78\right)$

其中，

Φ(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种类的机型中的某一个中发生了启动/关闭的变化。

因而，如下面那样。

[数83]

1～k_t的电设备的任一个变化为启动

1～k_t的电设备的任一个变化为关闭

无变化

问题在于估计该变化是k_t个电设备中的哪一个电设备的变化。为了进行该估计，如前所述，应用k-均值模型的考虑方法。其结果是，nzu_t ^k,k＝1,…,k_t，关于在给定u_t ⁰时它是基于哪个电设备的启动/关闭的变化而产生的后验概率。

那么，在给定u_t ⁰时，即使以nzu_t ^k,k＝1,…,k_t估计出关于k_t个电设备中的哪个电设备发生了变化的概率，也未必明确它给状态变量带来什么样的变化。

这是由于，根据发生变化之前的t-1期的状态变量的值x_t-1 ^k，状态变量的变化的方法不同。一般将根据t-1期的状态的不同而t期的状态转移的方法不同的模型称为马尔可夫转换(Markov-Switching)模型。

在本实施例中的用于NIALM的实时逆估计算法中，在状态转移中应用了该马尔可夫转换模型的考虑方法。

下面对此进行说明。

首先，从u_t ⁰>0的情况开始看。图5中示出了对u_t ⁰>0的情况下的状态变化进行总结而得到的表。

由于u_t ⁰>0，可认为k＝1,…,k_t的某一个设备存在启动的变化。但是，该变化根据t-1期的状态变量x_t-1 ^k的值而不同。

该变化一般来说能够如图5的表那样表达为转移概率的形式。即，将在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的机型、即已启动的设备未发生变化，为了表达这一点，在图5的表中，记入为P₁₁＝1、P₁₀＝0。

另一方面，对于在t-1期中x_t-1 ^k的值为0的机型而言，有可能发生启动的变化、即在t期中x_t ^k有可能变化为1，因此以nzu_t ^k来评价该概率P₀₁，这是本模型的构想。

按照该构想，如下那样设想转移概率。当设

[数84]

$\begin{array}{c}{p}_{01}=P(X_t^k=1|X_{t-1}^k=0)={\mathrm{nzu}}_t^k\\ p_{00}=1-{\mathrm{nzu}}_t^k\end{array}...\left(79\right)$

[数85]

$\begin{array}{c}{p}_{11}=P(X_t^k=1|X_{t-1}^k=1)=1\\ p_{10}=0\end{array}...\left(80\right)$

时，根据马尔可夫链(Markov chain)，x_t ^k取状态1的概率、即在t期中机型k的电设备变为启动的概率为下式。

[数86]

$P(X_t^k=1)=P(X_{t-1}^k=0)p_{01}+P(X_{t-1}^k=1)p_{11}...\left(81\right)$

现在，若着眼于随机变量的期望值表示概率，则下式(82)～(85)成立。

[数87]

$x_t^k=E\left(X_t^k\right)...\left(82\right)$

[数88]

$P(X_t^k=1)=E\left(X_t^k\right)=x_t^k...\left(83\right)$

[数89]

$P(X_{t-1}^k=1)=E\left(X_{t-1}^k\right)=x_{t-1}^k...\left(84\right)$

[数90]

$P(X_{t-1}^k=0)=1-x_{t-1}^k...\left(85\right)$

结果，当使用以上的关系时，式(82)成为状态方程，成为

[数91]

$x_t^k=(1-x_{t-1}^k){\mathrm{nzu}}_t^k+x_{t-1}^k+w_t^k...\left(86\right)$

[数92]

$x_t^k=x_{t-1}^k+(1-x_{t-1}^k){\mathrm{nzu}}_t^k+w_t^k...\left(87\right)$

，因此当再次将输入变量设为

[数93]

$u_t^k=(1-x_{t-1}^k){\mathrm{nzu}}_t^k...\left(88\right)$

时，能够如下那样表达为线性标准形

[数94]

$x_t^k=x_{t-1}^k+u_t^k+w_t^k...\left(89\right)$

同样地，在u_t ^k<0时为以下。在图6中示出了对u_t ⁰<0的情况下的状态变化进行总结而得到的表。

如图6的表的记载那样，如以下那样设想转移概率。

[数95]

$\begin{array}{c}{p}_{10}=P(X_t^k=0|X_{t-1}^k=1)={\mathrm{nzu}}_t^k\\ p_{11}=1-{\mathrm{nzu}}_t^k\end{array}...\left(90\right)$

[数96]

$\begin{array}{c}{p}_{01}=P(X_t^k=1|X_{t-1}^k=0)=0\\ p_{00}=1\end{array}...\left(91\right)$

根据马尔可夫链，x_t ^k取状态0的概率、即在t期中机型k的电设备变为关闭的概率如下。

[数97]

$P(X_t^k=0)=P(X_{t-1}^k=0)p_{00}+P(X_{t-1}^k=1)p_{10}...\left(92\right)$

另一方面，状态方程成为

[数98]

$1-x_t^k=(1-x_{t-1}^k)+x_{t-1}^k{\mathrm{nzu}}_t^k+w_t^k...\left(93\right)$

[数99]

$x_t^k=x_{t-1}^k-{\mathrm{nzu}}_t^k{x}_{t-1}^k+w_t^k...\left(94\right)$

，因此，当再次将输入变量u_t ^k设为

[数100]

$u_t^k=-{\mathrm{nzu}}_t^k{x}_{t-1}^k...\left(95\right)$

时，能够如下那样表达。

[数101]

$x_t^k=x_{t-1}^k+u_t^k+w_t^k...\left(96\right)$

另一方面，在u_t ⁰＝0时，

[数102]

$p(X_t^k=X_{t-1}^k)=p(u_t^k=0)=1...\left(97\right)$

即，转移概率为P₀₀＝1、P₁₁＝1。

因而，状态方程为下式。

[数103]

$x_t^k=x_{t-1}^k+w_t^k...\left(98\right)$

以上，能够导出卡尔曼滤波器的状态方程。在马尔可夫转换模型中，状态方程不是表达为单一的方程的形式，而是根据期的状态变量的状态、并同时根据总使用电力的跳跃为正、负或者0的情况，表达为多个方程的形式。

进一步附带地说，通常，在卡尔曼滤波器中，状态变量矢量的维数不发生变动而是固定的。但是，在本实施例的逆估计算法中，在时间经过的中途实时地估计个体电设备的机型和标准电力本身，机型数随着时间的经过而变化。因而，描述状态方程的状态变量矢量的维数根据时间的推移而变化，从而还出现状态变量矢量的维数变化时的变量的传递等未开拓的课题。

接着，使用以马尔可夫转换模型描述的状态方程和总使用电力的观测值，来说明利用卡尔曼滤波器的状态变量的最佳估计的递归性更新处理。

该处理是与图3的区域C-2相当的部分。本模型的状态变量是估计出的个体电设备的工作概率，因此得到总使用电力的观测值来随时更新个体电设备的工作状态的最佳估计。

状态变量的递归性的最佳更新算法是卡尔曼滤波器的核心。关于该计算处理，已针对卡尔曼滤波器的标准形描述为式(95)至式(65)。

卡尔曼滤波器是最小二乘法的在两个含义上的扩展。

第一点是，通常的最小二乘法将作为未知参数的β当作随机变量，求出其最小方差估计量，即，在初始时间点事先具有与其协方差有关的先验信息。

第二点是，在最小方差估计法中，在逐次输入了新的观测值时如何更新此前的β的估计值，即，将新的观测值所具有的、对向此前的估计值所延伸的空间的正交分量有贡献的部分加到现有的估计值来逐次进行更新。该第二点在卡尔曼滤波器中表达为使用被称为卡尔曼增益的项的状态变量的更新算法。

那么，试着依次将该计算处理变形为矩阵形式。

式的误差协方差R_t,t＝0,1,2,…是外生性地给定的。状态变量及其协方差矩阵x₀ ^{^+}、P₀ ⁺在t-1以后依次内生性地被修正，与此相对，Q_t、R_t在所有期中均是外生性地给定的。

接着，状态变量的右上标的区别是重要的。若再次示出它们，则为

[数104]

${\hat{x}}_t^{-}=E\left(x_t\right|y_1,y_2,...,y_{t-1})...\left(53\right)$

[数105]

${\hat{x}}_t^{+}=E\left(x_t\right|y_1,y_2,...,y_{t-1},y_t)...\left(54\right)$

(参照式(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))

[数106]

x_t＝F_t-1x_t-1+G_t-1u_t-1+w_t-1 …(48)

当将本逆估计算法具体地以矩阵形式写下去时，由于F_t-1、G_t-1为单位矩阵，因此如下。

[数107]

接着，按照状态方程的预测式(62)，来计算基于到t-1期为止的信息的、状态变量x_t的协方差P_t的估计值P_t ^-。

[数108]

$P_t^{-}=F_{t-1}{P}_{t-1}^{+}{F}_{t-1}^T+Q_{t-1}...\left(62\right)$

并且，根据下式来计算卡尔曼增益。

[数109]

$K_t=P_t^{-}{H}_t^T{(H_t{P}_t^{-}{H}_t^T+R_t)}^{-1}...\left(63\right)$

其中，H_t是由到t时间点为止估计出的个体电设备的标准电力构成的矢量，具体地说，如以下那样表达。

[数110]

$H_t=\left(\begin{array}{cccc}{P}_t^0& P_t^1& ...& P_t^{k_t}\end{array}\right)...\left(100\right)$

接着，观测t期的观测值、即总使用电力Pw_t，以t期的总使用电力的预测值H_tx_t ^{^-}与观测值之间的误差y_t-H_tx_t ^{^-}为基础，利用卡尔曼增益来修正状态变量x_t ^{^-}，由此如下那样求出使用了到t期为止的信息的、状态变量x_t的估计值x_t ^{^+}，其中，t期的总使用电力的预测值H_tx_t ^{^-}是基于通过使用直到t-1期为止的信息而得到的t期中的个体电设备的工作状态x_t ^{^-}与t期的个体电设备的标准电力P_t ^k的乘积和而得到的。

[数111]

${\hat{x}}_t^{+}={\hat{x}}_t^{-}+K_t(y_t-H_t{\hat{x}}_t^{-})...\left(101\right)$

如果将它表达为矩阵形式则为以下。

[数112]

$\left(\begin{array}{c}{\hat{x}}_t^{0+}\\ {\hat{x}}_t^{1+}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\hat{x}}_t^{k_t+}\end{array}\right)=\left(\begin{array}{c}{\hat{x}}_t^{0-}\\ {\hat{x}}_t^{1-}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\hat{x}}_t^{k_t-}\end{array}\right)+\left(\begin{array}{c}{K}_t^0\\ K_t^1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ K_t^{k_t}\end{array}\right)(y_t-\left(\begin{array}{cccc}{P}_t^0& P_t^1& ...& P_t^{k_t}\end{array}\right)\left(\begin{array}{c}{\hat{x}}_t^{0-}\\ {\hat{x}}_t^{1-}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\hat{x}}_t^{k_t-}\end{array}\right))...\left(102\right)$

最后，对于状态变量的估计量的最小协方差矩阵P的估计值，按下式那样从使用了到t-1期为止的信息的估计值P_t ^-更新为使用了到t期为止的信息的估计值P_t ⁺。

[数113]

$P_t^{+}=P_t^{-}-K_t{H}_t{P}_t^{-}...\left(103\right)$

通过以上的处理，作为期的输出，输出x_t ^{^+}、P_t ⁺以及H_t，并转移至t+1期。这是卡尔曼滤波器的第期的计算处理。

图7是将此前说明的本实施例的用于NIALM的实时逆估计算法中的整体的计算处理总结为流程图而得到的图。

在此，按照图7，在下面概述再次总结实时逆估计算法的流程而得到的内容。

1)给定初始值。在初始值中，给定个体电设备的工作概率x₀ ^{^+}及其协方差P₀ ⁺，还给定状态方程的误差协方差Q_t,t＝0,1,…以及观测方程的误差协方差R_t,t＝0,1,2,…。状态方程和观测方程的误差协方差在所有期中均是外生性地给定的。

2)通过k-均值算法，判断总使用电力的跳跃是否为此前存在的电设备的标准电力，如果是新的跳跃，则作为新的电设备追加到机型列表，将该标准电力临时登记为新机型的标准电力。

3)将临时登记的个体电设备的机型和标准电力包括在内地重新计算登记到机型列表中的电设备的标准电力及其方差。如果标准电力的重新计算的结果是在现有的机型之间、或者在新机型与现有的机型中有能够视作相同的，则将该电设备从机型列表去除。由此，将个体电设备的机型列表、其标准电力列表修正为最新的列表。

4)基于k-均值算法的计算结果，计算关于总使用电力的跳跃是最新的机型列表、标准电力列表中登记的、具有哪个标准电力的哪个个体电设备的工作状态的变化所引起的跳跃的发生难易度的概率。

5)构成根据总使用电力的跳跃是正是负还是0而状态转移不同的马尔可夫转换模型，构成根据总使用电力的跳跃的符号和t-1期的个体电设备的工作状态的不同而向t期的状态转移不同的状态方程。

6)利用使用了到t-1期为止的信息的个体电设备的工作概率以及t期中的个体电设备的标准电力的估计值的信息来预测t期的总使用电力。基于该预测值与t期的总使用电力的观测值的误差，通过卡尔曼滤波器来进行状态变量和状态变量的协方差矩阵的最佳估计的递归性更新，输出使用了到t期为止的信息的t期的个体电设备的工作概率x_t ^{^+}及其协方差矩阵P_t+。

本实施例的实时逆估计算法具有下面的特征和意义。

1)仅着眼于层使用电力的时间变动的1级级差、即跳跃。

2)进行以下简化：根据计划电力的不同来识别电设备的机型。

3)将跳跃作为针对簇的数量发生变化的可变的k-均值模型的输入，使用EM算法的解法，由此能够进行以1分钟为间隔的实时的对机型的识别和估计。

4)将可变的k-均值模型与状态空间建模的卡尔曼滤波器相结合，能够实时且动态地估计所识别出的个体电设备的工作状态的变化。

5)在k-均值模型和卡尔曼滤波器的合并中，以马尔可夫转换的考虑方法为强力的手段。

这样，本实施例的逆估计算法是意图能够进行实时的动态估计的实用性高的、到目前为止没有的算法。

接着，将本实施例的实时逆估计算法应用到实际的过程环境数据，来对该算法进行了评价。

如以下那样应用到实际的过程环境数据。

在实际生活的环境(东京都内)中将NIALM设置了8天(2013年3月12日～19日)，观测总使用电力，关于是否能够估计个体电设备的种类和标准电力进行了实验。NIALM使用模拟家庭环境以及与测量设备相同的设备，每隔1秒收集一次数据。协助实验的家庭是4人家庭，其中2人为小孩。实验期间的天气如图8所示的表。

在本实施例中，设想今后普及到一般家庭，使用了8天中的以1分钟为间隔的11520个数据。

在图9～图16中示出了该评价的结果。

图9示出了按总使用电力的1分钟的跳跃的规模的频数分布。根据该图明确可知，观测出的跳跃以0为中心在正和负对称地被描绘。因此，可以说，总使用电力的跳跃说明其设备变为启动、关闭的可能性非常高，该跳跃是标准电力的可能性高。

作为本实施例中构建的逆估计算法的应用对象，设想实际家庭环境中的总使用电力的以1分钟为间隔的时间变动数据。特别是，着眼于总使用电力的以1分钟为时间间隔的时间变动数据的1级级差、即跳跃，将其作为逆估计算法的输入流。

它是如下的考虑方法：在个体电设备被启动时，发生标准电力的正跳跃，在某一时间点个体电设备变为关闭时，发生负跳跃，个体电设备的使用电力变为0。该方法不像专利文献1所记载的技术那样还考虑总使用电力的波形，而是仅着眼于启动和关闭的跳跃的动作。

图10示出了按总使用电力规模的在总电力使用量中的占有率。当求出该值的合计时，得到总电力使用量8,941,000Wm。当将其换算为平均每个月的总电力使用量时，558.8kWh/月＝8941kWm÷60分钟÷8天×30天。可知，这接近平均每户的电力使用量544kWh/月(参照：(2013年3月)总务省家庭收支情况调查)的值，该被验者的家庭采取平均性的电使用方法，适合于NIALM的电力估计的实验对象。

图11是8天实际环境数据的实时处理算法的结果，是从开始使用起历经1200分钟的结果。在该图中横轴表示时间，纵轴表示总使用电力。

此外，在该算法中，未进行状态变量的传递。

图12是剪切出图11中的从开始使用起的150分钟并放大而得到的。另外，图13是剪切出图11中的150分钟到300分钟并放大而得到的。

根据图12可知，从开始使用起的150分钟内，存在几处在总使用电力的观测值与估计值之间产生误差的位置，但是在之后的150分钟以后，如图13那样，在观测值与估计值之间未产生大的误差，估计精度佳。

图14是使用以手动启动和关闭个体电设备的模拟家庭环境下的4小时内的以1分钟为间隔的数据并应用本实施例来估计个体电设备的种类和标准电力、并估计总使用电力和个体电设备的工作状态而得到的结果。在本实施例的卡尔曼滤波器中，将状态变量矢量的维数设为可变，但是在状态变量发生变化的情况下，需要考虑在变化后如何继承变化以前的状态变量。

图14未考虑个体电设备的种类的数量发生变化时的状态变量的传递，将个体电设备的工作概率清零。从开始使用起的150分钟内，个体电设备的种类、标准电力是未知的，尽管个体电设备的种类发生变化，但没有考虑状态变量的传递，因此存在好几个总使用电力的观测值与估计值大不相同的时间点。

但是，在150分以后，即使不考虑个体电设备的种类的数量发生变化时的状态变量的传递，个体电设备的种类的数量也通过学习而达到稳定的状态，总使用电力的观测值与估计值之间未产生大的差异。

与此相对，图15是估计个体电设备的种类和标准电力、并在考虑了个体电设备的种类的数量发生变化时的状态变量的传递的基础上估计总使用电力和个体电设备的工作状态而得到的结果。

与图13的未考虑状态变量的传递的情况相比，从开始使用起在总使用电力的观测值与估计值之间未出现大的差异，从初始的阶段起就实现了预测精度的提高。

如以上所说明的那样，实施例1的个体电设备工作状态估计装置及其方法具有以下的效果。

即，在实施例1的个体电设备工作状态估计装置及其方法中，即使不对各个电设备设置传感器、而且事先不具有存在几台什么种类的电设备之类的信息，也能够实时地监视各个电设备的使用状况(存在什么电设备的机型，哪个电设备处于什么样的工作状态)。因而，变得廉价，也容易应用于一般家庭等，其结果，能够实施电设备的使用状况的改善、电力供给控制的改善等来实现节能化。

另外，在由电设备机型确定部6a进行的聚类中使用k-均值算法，使在该k-均值算法中使用的与所述个体电设备的机型对应的簇的数量k可变，因此，即使中途追加新的电设备或者中途废弃此前使用的电设备从而作为监视对象的电设备的机型变得不同，也无需每次都进行机型列表的数据的追加输入、删除，从而也不需要其工时。其结果，即使是外行人也能够利用本实施例的个体电设备工作状态估计装置。

另外，k-均值算法使用具有期望值步骤和最大化步骤的EM算法，基于所输入的所述跳跃电力来估计所述跳跃电力为来自未知的k个正态分布的实现值时的k个正态分布的均值和方差，判定所输入的所述跳跃电力是来自哪个正态分布的实现值，因此，无论电设备增加或减少，即使不输入先验信息，也能够针对个体电设备获得所需的信息。

另外，EM算法中的期望值步骤和最大化步骤的每隔时间变动的时间间隔的迭代计算通常将迭代计算进行几次，但是在本实施例中仅进行一步重新计算，因此，根据EM算法的似然度的单调性，随着时间经过而估计精度提高，从而能够在充分维持估计精度的同时还将计算成本抑制得低。

另外，通常，在k-均值算法中使用全部数据，但是在本实施例的最大化步骤中，将使用的过去的数据的长度、返回的长度指定为比它之前得到的过去的数据少的长度，因此，不使用此前得到的全部数据就能够进行计算，从而能够将计算成本抑制得低。

另外，在EM算法中，在期望值步骤的期望值的计算中引入基于设该期望值的后验概率同所述期望值的先验分布与似然度之积成正比的贝叶斯学习的先验分布，来用于根据总使用电力和跳跃电力的过去的数据而得到的、个体电设备的启动或关闭的发生概率的时间变动、因使用环境的不同引起的变动、使用日变动中的至少一个变动的模式的学习结果的估计，因此能够进行最大似然估计。

另外，具备马尔可夫转换模型和卡尔曼滤波器，其中，马尔可夫转换模型与各电设备的工作概率的状态变量对应，卡尔曼滤波器使状态方程的标准形具有使所述跳跃电力与所述总使用电力对应的各输入变量，并使用作为所述各电设备的工作概率的状态变量，因此，能够使用马尔可夫转换模型对k-均值下的估计结果进行变换以使得能够在卡尔曼滤波器中利用该估计结果。其结果，能够动态且实时地推测个体电设备的工作状态。

另外，作为卡尔曼滤波器的状态变量的状态变量矢量的维数通常不变，但是使本实施例中的作为卡尔曼滤波器的状态变量的状态变量矢量的维数可变，因此，即使在中途追加新的电设备或者中途废弃此前使用的电设备的情况下，也无需每次都进行机型列表的数据的追加输入、删除，从而也不需要其工时。

以上，基于上述实施例说明了本发明，但是本发明不限于上述实施例，在不脱离本发明的宗旨的范围内进行的设计变更等也包含于本发明。

例如，本发明的个体电设备工作状态估计装置及其方法在上述实施例中应用于1户一般家庭，但是也可以不是一般家庭，或者也可以应用于将集体住宅集中而成的住宅。

或者，作为本发明的个体电设备工作状态估计装置及其方法的监视对象，也可以将1户房屋分为多个区域，对各个区域设置配电盘等，个别或综合地应用。

另外，电力计不限于安装于通常家庭的电力计，也可以利用其它电力计。另外，各处理中的运算不限于上述运算式。

另外，在上述实施例中，根据总使用电力的1级级差来求出跳跃电力，但是在发明中并不限于此。

产业上的可利用性

本发明的个体电设备工作状态估计装置及其方法能够利用于针对作为监视对象的多个电设备推测其工作状态的各种系统。

附图标记说明

1：房屋；2：布线；3a～3n：电设备；4：供给电线；5：电力计(总使用电力测定单元)；6：工作状态估计部；6a：时间变动计算部(时间变动计算单元)；6b：工作电设备机型确定部(工作电设备机型确定部)；6c：对应似然度估计部(对应似然度估计单元)；6d：个体电设备工作状态估计部(个体电设备工作状态估计单元)；6e：马尔可夫转换模型；6f：卡尔曼滤波器；7：显示器；8：配电盘(总使用电力测定单元)。

Claims (16)

1.一种个体电设备工作状态估计装置，其特征在于，具备：

总使用电力测定单元，其测定多个电设备的总使用电力；

时间变动计算单元，其基于由该总使用电力测定单元测定出的总使用电力来计算跳跃电力，该跳跃电力是该总使用电力的时间性变动；

电设备机型确定单元，其进行聚类，即，基于由该时间变动计算单元计算出的跳跃电力来判定是否要作为新的电设备追加到现有的机型列表，在要作为新的电设备追加到所述现有的机型列表的情况下，按个体电设备重新计算标准电力，根据需要而从所述现有的机型列表删除多余的机型，在不追加新的电设备的情况下维持所述现有的机型列表；

对应似然度估计单元，其基于由所述时间变动计算单元计算出的跳跃电力来估计关于升和降的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度；以及

个体电设备工作状态估计单元，其估计在输入了由该对应似然度估计单元得到的、在当前的个体电设备的工作状态下根据总使用电力的跳跃而得到的、与所述个体电设备的启动和关闭的事件发生有关的似然度时的所述个体电设备的工作状态的变化，对所述当前的个体电设备的工作状态的估计进行更新，根据使用该更新后的个体电设备的工作状态的估计结果和标准电力的估计结果所预测出的总使用电力的预测值以及所述总使用电力的测定数据，来优化所述个体电设备的机型的估计、所述标准电力的估计以及所述工作状态的估计，以测定出的所述总使用电力的时间变动的时间间隔来动态地估计时刻变化的所述个体电设备的工作概率。

2.根据权利要求1所述的个体电设备工作状态估计装置，其特征在于，

在由所述电设备机型确定单元进行的聚类中使用k-均值算法，使在该k-均值算法中使用的与所述个体电设备的机型对应的簇的数量k可变。

3.根据权利要求2所述的个体电设备工作状态估计装置，其特征在于，

所述k-均值算法使用具有期望值步骤和最大化步骤的EM算法，基于所输入的所述跳跃电力来估计所述跳跃电力为来自未知的k个正态分布的实现值时的k个正态分布的均值和方差，判定所输入的所述跳跃电力是来自哪个正态分布的实现值。

4.根据权利要求3所述的个体电设备工作状态估计装置，其特征在于，

所述EM算法中的期望值步骤和最大化步骤的每隔时间变动间隔的迭代计算仅进行一步重新计算。

5.根据权利要求3或4所述的个体电设备工作状态估计装置，其特征在于，

所述最大化步骤中使用的均值和方差的过去的数据的长度、返回的长度被指定为比该过去的数据之前得到的过去的数据少的长度。

6.根据权利要求3至5中的任一项所述的个体电设备工作状态估计装置，其特征在于，

所述EM算法在所述期望值步骤的期望值的计算中引入基于设该期望值的后验概率同所述期望值的先验分布与似然度之积成正比的贝叶斯学习的所述先验分布，来用于根据所述总使用电力和所述跳跃电力的过去的数据而得到的、所述个体电设备的启动或关闭的发生概率的时间变动、因使用环境的不同引起的变动、使用日变动中的至少一个变动的模式的学习结果的估计。

7.根据权利要求2至6中的任一项所述的个体电设备工作状态估计装置，其特征在于，具备：

马尔可夫转换模型，其与各电设备的工作概率的状态变量对应；以及

卡尔曼滤波器，其使状态方程的标准形具有使所述跳跃电量与所述总使用电力对应的各输入变量，并使用作为所述各电设备的工作概率的状态变量。

8.根据权利要求7所述的个体电设备工作状态估计装置，其特征在于，

使所述卡尔曼滤波器的状态变量为维数可变的状态变量矢量。

9.一种个体电设备工作状态估计方法，其特征在于，具有以下步骤：

测定多个电设备的总使用电力；

基于所测定出的所述总使用电力来计算跳跃电力，该跳跃电力是该总使用电力的时间性变动；

进行聚类，即，基于所计算出的所述跳跃电力来判定是否要作为新的电设备追加到现有的机型列表，在要作为新的电设备追加到所述现有的机型列表的情况下，按个体电设备重新计算标准电力，根据需要而从所述现有的机型列表删除多余的机型，在不追加新的电设备的情况下维持所述现有的机型列表；

基于所计算出的所述跳跃电力来估计关于升和降的跳跃与哪个电设备的启动和关闭的事件的发生相对应的似然度；以及

估计在输入了在当前的个体电设备的工作状态下根据总使用电力的跳跃而得到的、与所述个体电设备的启动和关闭的事件发生有关的似然度时的该个体电设备的工作状态的变化，对所述当前的个体电设备的工作状态的估计进行更新，根据使用该更新后的个体电设备的工作状态的估计结果和标准电力的估计结果所预测出的总使用电力的预测值以及所述总使用电力的测定数据，来优化所述个体电设备的机型的估计、所述标准电力的估计以及所述工作状态的估计，以测定出的所述总使用电力的时间变动的时间间隔来动态地估计时刻变化的所述个体电设备的工作概率。

10.根据权利要求9所述的个体电设备工作状态估计方法，其特征在于，

所述聚类使用k-均值算法，使在该k-均值算法中使用的与所述个体电设备的机型对应的簇的数量k可变。

11.根据权利要求10所述的个体电设备工作状态估计方法，其特征在于，

所述k-均值算法使用具有期望值步骤和最大化步骤的EM算法，基于所输入的所述跳跃电力来估计所述跳跃电力为来自未知的k个正态分布的实现值时的k个正态分布的均值和方差，判定所输入的所述跳跃电力是来自哪个正态分布的实现值。

12.根据权利要求11所述的个体电设备工作状态估计方法，其特征在于，

所述EM算法中的期望值步骤和最大化步骤的每隔时间变动间隔的迭代计算仅进行一步重新计算。

13.根据权利要求11或12所述的个体电设备工作状态估计方法，其特征在于，

所述最大化步骤中使用的均值和方差的过去的数据的长度、返回的长度被指定为比该过去的数据之前得到的过去的数据少的长度。

14.根据权利要求11至13中的任一项所述的个体电设备工作状态估计方法，其特征在于，

所述EM算法在所述期望值步骤的期望值的计算中引入基于设该期望值的后验概率同所述期望值的先验分布与似然度之积成正比的贝叶斯学习的所述先验分布，来用于根据所述总使用电力和所述跳跃电力的过去的数据而得到的、所述个体电设备的启动或关闭的发生概率的时间变动、因使用环境的不同引起的变动、使用日变动中的至少一个变动的模式的学习结果的估计。

15.根据权利要求9至14中的任一项所述的个体电设备工作状态估计方法，其特征在于，具有以下步骤：

使用马尔可夫转换模型来与各电设备的工作概率的状态变量对应；以及

使卡尔曼滤波器的状态方程的标准形具有使所述跳跃电量与所述总使用电力对应的各输入变量，在所述卡尔曼滤波器中使用作为所述各电设备的工作概率的状态变量。

16.根据权利要求15所述的个体电设备工作状态估计方法，其特征在于，

使所述卡尔曼滤波器的状态变量为维数可变的状态变量矢量。