CN104505894B - 一种基于矿用锂离子电池的电源管理系统及状态估计方法 - Google Patents

Abstract

一种基于矿用锂离子电池的电源管理系统及状态估计方法。包括两部分：(1)基于多节单体锂电池组成的矿用锂离子电池组的电源管理系统，包括主控单元、单总线温度采集单元、电源管理单元、电流采集单元、控制单元和数据传输单元，实现过欠压保护、过流保护、超温保护、充电均衡管理、容量自检及运行参数的实时上传和修改功能；(2)基于高斯过程推理学习方法的状态估计方法，把稀疏高斯过程和高斯过程假定密度滤波器，考虑了策略噪声、系统噪声和模型的不确定性，使用解析的方法以较少的计算代价进行滤波、预测运算，易于对非平稳动态系统建模。应用效果表明，电源管理系统和状态估计模型有效地保证矿用电池安全高效工作，最大化电池能量利用率。

Description

一种基于矿用锂离子电池的电源管理系统及状态估计方法

技术领域

本发明涉及一种电源管理系统及状态估计方法，特别是一种基于矿用锂离子电池的电源管理系统及状态估计方法。

背景技术

随着煤炭工业发展和矿山装备技术进步，监测通信系统、紧急避险设施、井下辅助运输设备等煤矿装备对储能电源的要求越来越高，特别是灾变环境下，由于安全原因一般都会失去外电供应，依靠电池供电的备用电源成为首选。而锂离子电池凭借其以其体积小、能量高、输出电流大、效率高、高安全可靠性、无污染等成为矿用电池的首选。

目前，对由多个锂离子电池单体串联构成的矿用电源进行电源管理，实际应用在矿用监控设备像信号采集装置等的供电电源的电源管理系统几乎没有或不完善，使得供电电源供电电压不稳定，容量损耗大，过充、过放等，这严重影响安全性和使用寿命，电池组充电频率高、更换电池任务繁重等也给工人增加了工作强度。

在实际使用时，电池组容量会随着循环次数衰减，同时环境温度、充放电倍率也会显著影响充电状态(State of Charge，SOC)。一般而言，电池放电深度也会严重影响电池的容量和寿命，锂电阻变化与电池寿命影响密切。在充电过程中，电流和电压会随着时间的变化而变化，能够很好表征充电模式。而充电状态/容量与上述因素是复杂的动态非线性关系。为了建立复杂动态非线性系统以估计电池容量和电池剩余寿命，人们做了如下几个方面研究：

(1)电池充电、放电、寿命的机理模型研究。该方面包括电池充电过程中电压、电流、以及电池部件的变化关系，不同因素如搁置时间、循环次数、温度、阻抗、充放电倍率、放电深度、电化学特征等因素对电池容量衰退趋势影响，结合电池动力学特征，从而建立等价电路模型以及多种数学模型。该类方法具有形式简单、工作机理明确、求解容易等优点。然而该方法无法描述电池工作过程中出现的物理和化学变化，因而对SOC、SOH等电池重要参数无法预测。

(2)动态非线性动态模型的辨识。常见的辨识方法有Kalman滤波、扩展Kalman滤波方法、粒子滤波方法等非线性动态系统建模方法。基于Kalman滤波方法虽然具有明确的解析形式，但是该方法只适用的线性状态空间建模。基于粒子滤波的非线性状态空间建模方法需要付出很大的代价，并且算法的性能受到粒子数量、粒子退化等因素影响。

(3)模糊理论的建模方法。根据电池电极电化学测量输出或者根据不同工作温度、不同健康状态条件下的电池内阻变化建立模糊逻辑模型，以此辨识等价的电路电阻模型。还有其它一些电池状态估计方法，比如Bayesian Mente Carlo方法自回归积分滑动均值(ARIMA)、相关支持向量机(RVM)、支持向量机(SVM)等电池状态预测方法。

发明内容

本发明的目的是要提供一种能有效地保证矿用电池安全高效工作，最大化电池能量利用率的基于矿用锂离子电池的电源管理系统及状态估计方法。

本发明的目的是这样实现的：本发明包括电池电源管理系统和电池状态估计方法：

所述的电池电源管理系统包括电源管理装置和电源管理方法，所述的电源管理装置是基于多节单体锂电池组成的矿用锂离子电池组进行管理的电源管理系统，包括主控单元、单总线温度采集单元、电源管理单元、电流采集单元、控制单元和数据传输单元；主控单元的输入端通过数据线与单总线温度采集单元、电流采集单元和数据传输单元的输出端连接，主控单元的输出端通过数据线与电源管理单元和控制单元的输入端连接。

所述的主控单元采用增强型系列微控制器STM32F103R8T6构成。

所述的单总线温度采集单元基于单线数据传输的温度传感器DS18B20构成。

所述的电源管理单元，基于锂离子电池监控芯片AD7280A构成。

所述的电源管理方法：主控单元和单总线温度采集单元，在电源管理系统首次上电后，主控单元加电自检并配置电源管理单元的参数，参数配置完成后，系统等待具有唯一序列号的单总线温度传感器接入，主控单元间隔搜索单总线上的温度传感器数量及各自的序列号，并存储在具有掉电不丢失数据特性的FLASH中，此处是由主控单元的FLASH模拟EEPROM实现；

所述的电流采集单元，根据采集的数据和系统以及不同电流等级和精度要求，共分为两组不同的电流采集方案：a、充电电流采集采用霍尔电流采集隔离方案；b、放电和容量自检电流采样通过串联精密电阻和电压变换测量，电流采集是基于主控单元自带的12位AD模块，外部加载独立高精度的电压源，结合特定的采样方案和数据处理方法，实现高精度的电流采样；

所述的电源管理单元，主控单元通过SPI接口操作电源管理单元，实现电压采集和充电均衡管理任务，经任务调度，可以不停的回读单体电池的电压、电池组的总电压；充电时，根据单体电池不同的电压等级，主控单元可以选择不同的均衡阀值来进行充电均衡管理，同时主控单元判断各部分的电压采集是否正常，以此判断电压采集线是否出问题，并更新温度和状态数据；

所述的数据传输单元，通过RS232接口与主控单元进行通信，接收电池数据，用于数据的显示，并给主控单元发出参数设置指令，可以设置调整的参数有：单体电池过压阀值、单体电池欠压阀值、充放电电流阀值、充电均衡时间、过流保护恢复时间、充电均衡电压、电池组容量值、充放电和容量自检电流修正系数。

所述的电池状态估计方法是基于高斯过程推理学习方法的状态估计方法，该方法是非线性状态空间建模方法，把稀疏高斯过程和高斯过程假定密度滤波器，既考虑了策略噪声，又考虑了系统噪声和模型的不确定性，使用解析的方法以较少的计算代价进行滤波、预测运算，易于对非平稳动态系统建模。

有益效果，由于采用了上述方案，利用基于多节单体锂电池组成的矿用锂离子电池组的电源管理系统，一方面实现了过欠压保护、过流保护、超温保护、充电均衡管理、容量自检，另一方面可以获取对充电过程、放电过程产生的非线性运行参数，通过采用基于高斯过程推理学习的状态估计算法可以对锂电池的运行状态进行有效的分类和准确的诊断。

本发明设计合理的电源管理系统，能够对电池外特性如电压、电流、温度等参数进行准确的数据采集，电池充放电均衡控制，电池剩余电量状态(SOC)的预测，电池过充、过放、过压等安全管理，并且有电池热管理等。这样可以有效地保证电池安全高效工作，避免过热、欠压等，最大化电池能量利用率。同时矿用电源经电源管理系统提供本质安全电源供电，使得本质安全型矿用设备使用很方便。

应用效果表明，电源管理系统和状态估计模型有效地保证矿用电池安全高效工作，最大化电池能量利用率，有效地延长电池使用寿命。

附图说明：

图1是本发明电源管理系统原理框图。

图2是本发明电源管理系统方案框图。

具体实施方式

本发明包括电池电源管理系统和电池状态估计方法：

所述的电池电源管理系统包括电源管理装置和电源管理方法，所述的电源管理装置是基于多节单体锂电池组成的矿用锂离子电池组进行管理的电源管理系统，包括主控单元、单总线温度采集单元、电源管理单元、电流采集单元、控制单元和数据传输单元；主控单元的输入端通过数据线与单总线温度采集单元、电流采集单元和数据传输单元的输出端连接，主控单元的输出端通过数据线与电源管理单元和控制单元的输入端连接；实现过欠压保护、过流保护、超温保护、充电均衡管理、容量自检以及运行参数的实时上传和修改的功能。

所述的主控单元采用增强型系列微控制器STM32F103R8T6构成；

所述的单总线温度采集单元基于单线数据传输的温度传感器DS18B20构成。

所述的电源管理单元，基于锂离子电池监控芯片AD7280A构成。

所述的电源管理方法：主控单元和单总线温度采集单元，在电源管理系统首次上电后，主控单元加电自检并配置电源管理单元的参数，参数配置完成后，系统等待具有唯一序列号的单总线温度传感器接入，主控单元间隔搜索单总线上的温度传感器数量及各自的序列号，并存储在具有掉电不丢失数据特性的FLASH中，此处是由主控单元的FLASH模拟EEPROM实现。

所述的电流采集单元，根据采集的数据和系统以及不同电流等级和精度要求，共分为两组不同的电流采集方案：a、充电电流采集采用霍尔电流采集隔离方案，霍尔电流传感器选用低偏移线性霍尔传感器专用芯片ACSTM12-20A；b、放电和容量自检电流采样通过串联精密电阻和电压变换测量，电流采集是基于主控单元自带的12位AD模块，外部加载独立高精度的电压源，结合特定的采样方案和数据处理方法，实现高精度的电流采样。

所述的电源管理单元，基于锂离子电池监控芯片AD7280A设计，主控单元通过SPI接口操作电源管理单元，实现电压采集和充电均衡管理任务，经任务调度，可以不停的回读单体电池的电压、电池组的总电压；充电时，根据单体电池不同的电压等级，主控单元可以选择不同的均衡阀值来进行充电均衡管理，同时主控单元判断各部分的电压采集是否正常，以此判断电压采集线是否出问题，并更新温度和状态数据。

所述的数据传输单元，通过RS232接口与主控单元进行通信，接收电池数据，用于数据的显示，并给主控单元发出参数设置指令，可以设置调整的参数有：单体电池过压阀值、单体电池欠压阀值、充放电电流阀值、充电均衡时间、过流保护恢复时间、充电均衡电压、电池组容量值、充放电和容量自检电流修正系数。

所述的电池状态估计方法是基于高斯过程推理学习方法的状态估计方法，该方法是非线性状态空间建模方法，把稀疏高斯过程和高斯过程假定密度滤波器，既考虑了策略噪声，又考虑了系统噪声和模型的不确定性，使用解析的方法以较少的计算代价进行滤波、预测运算，易于对非平稳动态系统建模。

实施例1：

下面结合附图对本发明的实施作进一步的说明：

本发明提供的电源管理系统是基于多节单体锂电池组成的电池组进行管理，电源管理系统原理框图如图1所示，采用如图2所示的技术方案。细节技术方案实现如下：

系统首次上电后，主控单元加电自检并配置电源管理单元的参数。

参数配置完成后，系统等待具有唯一序列号的单总线温度传感器接入。主控单元间隔搜索单总线上的温度传感器数量及各自的序列号，并存储进缓存。当总线上的传感器数量比现有的数量多1的时候，主控单元将前后两次搜索到的序列号对比，并提取出多出来的一组传感器序列号值。由此，主控单元可根据传感器接入系统的顺序区分不同的传感器。当传感器数量达到5组，即温度传感器全部接入系统。系统按照传感器接入总线的顺序整理相应的序列号，并记录在具有掉电不丢失数据特性的FLASH中，此处是由主控单元的FLASH模拟EEPROM实现。

系统进入主程序后，不同的任务由RTOS调度。同一时刻，只能运行对外输出、充电、和容量自检三个功能之一。不同的任务实现不同的功能，细节如下：

电压采集和充电均衡任务：该部分是由主控单元通过SPI接口操作电源管理芯片实现，经任务调度，可以不停的回读单体电池的电压、电池组的总电压。充电时，根据单体电池不同的电压等级，主控单元可以选择不同的均衡阀值来进行充电均衡。同时主控单元判断各部分的电压采集是否正常，以此判断电压采集线是否出问题，并更新温度和状态数据。

电流采集任务：根据采集到的数据和系统根据不同的电流等级和精度要求，共分为三组不同的电流采集方案。如图2所示，充电电流采集采用霍尔电流采集隔离方案，放电和容量自检电流采样通过串联精密电阻和电压变换测量。整个系统的电流采集是基于主控单元自带的12位AD模块，外部加载独立高精度的电压源，结合特定的采样方案和数据处理方法，实现高精度的电流采样。AD模块不仅采集三组电流值，同时采集基准端的模拟地和参考电压。按以下计算公式计算出三组电流采样对应的实际电压。

实际电压＝(实际采样值-模拟地采样值)*(参考电压-0V)/(参考电压采样值-模拟地采样值)进而根据电流电压变换关系，换算出实际电流。

温度采集任务：主控单元每隔20ms发送一次温度查询指令，循环读取各个单体电池的温度，同时系统判断挂接在总线上的温度传感器是否有问题，并更新温度和状态数据。

系统监测任务：该任务比较判断由电压电流采集任务和温度采集任务得来的数据，并做相应的处理。本安放电电流最大900mA，发生过流，主控单元将在80ms内切断与充电过流、放电过流和容量自检过流相对应的输出。发生输出短路，主控单元将在40ms内切断对外输出。温度保护阀值是60摄氏度，达到阀值后，将关断容量自检输出、充电和放电。达到55摄氏度，进行预报警。电压保护分为过压和欠压。系统达到欠压阀值，将关闭对外输出，同时打开充电端口。系统达到过压阀值，关闭充电端口，并打开对外输出和容量自检功能。主控单元检测到对外负载，关闭除此之外的其它端口。同时，主控单元还进行过充过放失效保护检测，监测电源管理系统的充放电控制单元是否失效。关断输出后，单体电池电压依然下降，说明输出关断失效。关断充电后，单体电池电压依然升高，说明充电控制端口失效。

数据传送和接收任务：主控单元通过图2中的RS232接口实现对外通信。实现主控单元对外传输电池数据，用于数据的显示，并根据接收到的指令实现参数设置。可以设置调整的参数有：单体电池过压阀值、单体电池欠压阀值、充放电电流阀值、充电均衡时间、过流保护恢复时间、充电均衡电压、电池组容量值、充放电和容量自检电流修正系数。

下面对电池状态估计方法做进一步说明：

1.稀疏高斯过程逼近

1.1高斯过程

给出随机变量集X＝{f(x)|x∈R^D}，高斯过程(GP)就是假定X上的任意子集均服从相同的高斯分布。作为一个鲁棒非参数Bayesian模型，GP只需由均值m(x)和协方差函数k(x,x')确定。假定给出随机函数f(x)的先验，关于GP输出的后验输出可以逼近任意的未知非线性函数，同时考虑了不确定输入。

考虑一个训练数据集D＝{(x_t,y_t),i＝1,2,...,N}，其中x_t为输入y_i为含噪标量输出，每个标量输出y_i可以由高斯噪声ε_i和在位置x_i处隐藏函数f(x_i)之和，即

y_i＝f(x_i)+ε_i (1)

这里ε_i～N(0,σ²)。记X为所有的训练输入组成的矩阵，y为所有输出组成的向量以及f为所有隐藏函数组成的向量。关于隐藏函数的先验为

p(f|x₁,x₂,...,x_N)＝N(0,K) (2)

这里K是N×N协方差矩阵，其第(i,j)元素为K_ij＝k(x_i,x_j)。广泛应用的协方差函数是一个平稳均方指数函数。

$k_{\mathrm{se}}(x_i,x_j)={\mathrm{\σ}}_{\mathrm{se}}^2\exp (-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{l}_d^2(x_{\mathrm{id}}-x_{\mathrm{jd}}))---\left(3\right)$

这里是尺度参数，l_d是长尺度用于控制d维方向的相关性下降速率。y和f联合分布为p(y,f)＝p(y|f)p(f)。根据Bayesian规则，后验均值函数m_y(x)和协方差函数k_y(x,x')由后验GP得出：

m_y(x)＝E_h[h(x*)]＝k_xN(σ²I+K_N×N)^-1y (4)

k_y(x,x')＝var_h[h(x*)]＝k(x,x')-k_xN(σ²I+K_N×N)^-1k_Nx (5)

这里K_N×N是关于训练输入的N×N协方差矩阵，k_xN＝[k(x,x₁),k(x,x₂),...,k(x,x_N)]为协方差向量并且这样可以计算关于训练隐藏函数的高斯后验p(f|y)。在输入点x*处的预测输出y*＝f(x*)+ε*可以写为

p(y*|y,X,x*)＝N(y*|m_y(x*),k_y(x*,x*)+σ²) (6)

这里 $k_y(x*,x*)=k(x*,x*)-k_{*,N}{K}_{N\×N}^{-1}{k}_{N,*}^{-1},k_{*,N}=[k(x*,x_1),k(x*,x_2),...,k(x*,x_N)],m_y(x*)=k_{*,N}{K}_{N\×N}^{-1}y.$

GP的最优模型参数可以通过最大化下面边缘似然函数对数求取，即

log p(y)＝log(N(y|0,σ²I+K_N×N)) (7)

GP需要计算N×N矩阵的逆，这样对大规模数据集需要花费很大的计算代价。为了解决这个问题，稀疏GP(SGP)通过逼近(7)所示真实似然函数推断超参数。SGP另外一个重要优点是能够出处理不确定输入变量，这个性质对使用不直接观测隐藏状态训练非线性转换函数非常关键。

1.2稀疏高斯过程

SGP对GP引入了条件独立性假设。这里简要介绍全独立训练条件(FITC)和部分独立训练训练条件(FITC)SGP。

稀疏高斯过程逼近是基于含有M(M＝N)个伪输入(也称为诱导输入)的小规模数据集和相关的伪输出u＝[u₁,u₂,L,u_M]^T，其Gaussian先验为

$u|\stackrel{\‾}{X}:N(0,K_{\mathrm{uu}})---\left(8\right)$

这里是协方差矩阵。训练观测集合隐含目标函数f，相应的似然函数为

$p\left(f\right|X,u,\stackrel{\‾}{X})=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(f\left(x_i\right)\right|x_i,u,\stackrel{\‾}{X})=N(K_{\mathrm{fu}}{K}_{\mathrm{uu}}^{-1}u,\mathrm{diag}(K_{\mathrm{ff}}-Q_{\mathrm{ff}}))---\left(9\right)$

这里 $Q_{\mathrm{ff}}=K_{\mathrm{fu}}{K}_{\mathrm{uu}}^{-1}{K}_{\mathrm{uf}},K_{f,u}={\left[k(x_i,{\stackrel{\‾}{x}}_j)\right]}_{i,j=1}^{N,M},K_{f,u}=K_{u,f}^T,K_{f,f}={\left[k(x_i,x_j)\right]}_{i,j=\mathrm{1,2}}^N.$ 因此有

$p\left(y\right|X,u,\stackrel{\‾}{X})=N\left(y\right|K_{f,u}{K}_{\mathrm{uu}}^{-1}u,\mathrm{diag}[K_{f,f}-Q_{f,f}]+{\mathrm{\σ}}_{\mathrm{noise}}^{2}I).$

FITC逼近方法假设训练集和测试变量完全条件独立，即

$p(f,f_{*}|X_{*},X,u,\stackrel{\‾}{X})=p\left(f_{*}\right|X_{*},u,\stackrel{\‾}{X})p\left(f\right|X,u,\stackrel{\‾}{X})X,X_{*}{\∈}^{N\×D}---\left(10\right)$

$p\left(f\right|X,u,\stackrel{\‾}{X})=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(f\left(x_i\right)\right|x_i,u,\stackrel{\‾}{X})=N(K_{\mathrm{fu}}{K}_{\mathrm{uu}}^{-1}u,\mathrm{diag}(K_{\mathrm{ff}}-Q_{\mathrm{ff}}))---\left(11\right)$

$p\left(f_{*}\right|X_{*},u,\stackrel{\‾}{X})=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(f\left(x_{*i}\right)\right|x_{*i},u,\stackrel{\‾}{X})=N(K_{f_{*}u}{K}_{\mathrm{uu}}^{-1}u,\mathrm{diag}(K_{f_{*}{f}_{*}}-Q_{f_{*}{f}_{*}}))---\left(12\right)$

这里 $Q_{\mathrm{ff}}=K_{\mathrm{fu}}{K}_{\mathrm{uu}}^{-1}{K}_{\mathrm{uf}},{Q_{f_{*}{f}_{*}}=K_{f_{*}u}{K}_{\mathrm{uu}}^{-1}{K}_{u{f}_{*}},K}_{f,u}={\left[k(x_i,{\stackrel{\‾}{x}}_j)\right]}_{i,j=1}^{N,M},K_{f,u}=K_{u,f}^T,K_{f,f}={\left[k(x_i,x_j)\right]}_{i,j=\mathrm{1,2}}^N.$ 使用u的先验对u进行积分，可以得到关于f和f_*的先验。

p(f|X)＝N(f|0,Q_ff+diag(K_ff-Q_ff)) (13)

$p\left(f_{*}\right|X_{*})=N\left(f_{*}\right|0,Q_{f_{*}{f}_{*}}+\mathrm{diag}(K_{f_{*}{f}_{*}}-Q_{f_{*}{f}_{*}}))---\left(14\right)$

令Σ_ff＝Q_ff+diag(K_ff-Q_ff)。FITC关于f和f_*的先验联合分布由下式可得

$p(f,f_{*}|X,X_{*})=N(0,{\mathrm{\Σ}}_{f,f_{*}})=N(0,\left[\begin{array}{cc}{\mathrm{\Σ}}_{\mathrm{ff}}& Q_{{\mathrm{ff}}_{*}}\\ Q_{f_{*}f}& K_{f_{*}{f}_{*}}\end{array}\right])---\left(15\right)$

从(15)和(9)可以得到基于训练输入X和诱导输入的边缘似然函数

$q\left(y\right|\stackrel{\‾}{X},X)=\∫p\left(y\right|f)p\left(f\right|X,u,\stackrel{\‾}{X})p\left(u\right|\stackrel{\‾}{X})\mathrm{dfdu}=N\left(y\right|0,K_{\mathrm{fu}}{K}_{\mathrm{uu}}^{-1}{K}_{\mathrm{uf}}+\mathrm{\Λ})---\left(16\right)$

这里Λ＝diag(K_ff-Q_ff)+σ²I。通过最大化下面对数似然函数求取最优伪输入位置和超参数

$\underset{\mathrm{\θ},\stackrel{\‾}{X},{\mathrm{\σ}}^2}{\arg \max }\log q\left(y\right|\stackrel{\‾}{X})=-\frac{1}{2}\log \left|\right|Q_{f,f}+\mathrm{\Λ}\left|\right|-\frac{1}{2}{y}^T{(Q_{f,f}+\mathrm{\Λ})}^{-1}y-\frac{N}{2}\log 2\mathrm{\π}---\left(17\right)$

使用关于和超参数θ的梯度实现上述似然函数优化。给出新的输入x*，与常规GP一样，从(6)同样可计算预测分布

$\begin{array}{c}p(y*|x*,\stackrel{\‾}{X},y,X)=\∫p(y*|x*,u,\stackrel{\‾}{X})p\left(u\right|y,X,\stackrel{\‾}{X})d\stackrel{\‾}{f}\\ =N(k_{*u}^T{\mathrm{\Σ}}^{-1}{K}_{u,f}{(\mathrm{\Λ}+{\mathrm{\σ}}^{2}I)}^{-1}y,k_{*,*}-Q_{*,*}+k_{*f}^T(K_{\mathrm{uu}}^{-1}-{\mathrm{\Σ}}^{-1})k_{f*}+{\mathrm{\σ}}^2)\end{array}---\left(18\right)$

这里Σ＝K_uu+K_u,f(Λ+σ²I)^-1K_f,u。

f(x_i)与f(x_j)的协方差可以写为

k_f(x_i,x_j)＝cov(f(x_i),f(x_j))＝E_f((f(x_i)-m_f(x_i))(f(x_j)-m_f(x_i))^T)

这里E_f表示对函数f(x)的期望，k_f(x_i,x_j)称为核函数。常用的协方差核函数是均方指数(SE)协方差函数，该函数包括自相关确定(ARD)和方差协方差，其形式为

$\begin{array}{c}{k}_f(x_i,x_j)=k_{\mathrm{se}}(x_i,x_j)+k_{\mathrm{noise}}(x_i,x_j)\\ =a^2\mathrm{ecp}(-\frac{1}{2}{(x_i-x_j)}^T\mathrm{\Λ}(x_i-x_j))+{\mathrm{\δ}}_{a,b}{\mathrm{\σ}}_{\mathrm{\ϵ}}^2\end{array}---\left(19\right)$

这里是尺度长矩阵，a²是函数f的方差。

2.高斯过程推理与学习

一般来说，非线性状态空间模型通常用于描述非线性动态系统。状态空间模型假定存在隐藏状态序列x_t，该隐藏序列根据转换函数(或系统函数)f:R^D→R^D确定的Markov过程随着时间不断演化。通过由y_t通过测量函数g:R^D→R^d间接观测隐藏函数。非线性状态方程模型的一般形式由下式给出

x_t＝f(x_t-1)+w_t,x_t∈R^D (20)

y_t＝g(x_t)+v_t,y_t∈R^d (21)

这里w_t～N(0,Σ_w)andv_t～N(0,Σ_v)分别是独立同分布的高斯系统噪声和测量噪声。离散时间步t从0到T。时间序列的初始状态x₀的分布是Gaussian先验分布实际上，高维时间序列可以由低维隐藏状态来描述。函数f和g都是未知的且由GP、建模，即f～GP_f，g～GP_g。需要预测步(从x_t-1移动到x_t)和滤波步(从观测y_t到隐藏状态x_t)逼近真实的隐藏状态

1)预测步

该步使用以前的滤波结果p(x_t-1|y_1:t-1)作为先验确定隐藏状态x_t的分布p(x_t|y_1:t-1)，这里下标1:t-1表示1,2,...,t-1的简写。根据Bayesian定理有

p(x_t|y_1:t-1)＝∫p(x_t-1|y_1:t-1)p(x_t|x_1:t-1)dx_t-1 (22)

上式无法解析地求解，因而需要近似方法。

2)滤波步骤

当f和g是未知的，并分别由带有SE核k_f、k_g的GPGP_f和GP_g建模时，GP-ADF能够计算近似高斯测量预测分布 $p\left(y_t\right|y_{1:t-1})\≈N({\mathrm{\μ}}_{t|t-1}^y,{\mathrm{\Σ}}_{t|t-1}^y),$ 状态预测分布 $p\left(x_t\right|y_{1:t-1})\≈N({\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)$ 以及带有不确定输入x_t的联合分布分布本专利使用简单标记当a＝μ表示均值，a＝Σ表示协方差，b表示考虑的时间步，c表示当前考虑测量的时间步，d∈{x,y}表示隐藏状态或者测量。

2.1测量预测

测量预测分布由如下高斯分布逼近

$p\left(y_t\right|y_{1:t-1})=\∫p(y_t,x_t|y_{1:t-1}){\mathrm{dx}}_t=\∫p\left(y_t\right|x_t)p\left(x_t\right|y_{1:t-1}){\mathrm{dx}}_t\≈N({\mathrm{\μ}}_{t|t-1}^y,{\mathrm{\Σ}}_{t|t-1}^y)---\left(23\right)$

使用式(21)测量模型和状态预测分布从y_1:t-1预测y_t。使用FITCSGP回归模型建立式(21)所示测量模型，状态空间模型由所学测量模型超参数伪输入ξ∈R^M×D以及输出υ∈R^M×E组成。假设给出不确定输入x_t|t-1：p(x_t|y_1:t-1)＝N(μ_t|t-1,Σ_t|t-1)，式(21)的每一维输出关于伪训练数据是条件独立的，即

g_a(x_t|t-1)⊥g_b(x_t|t-1)|x_t|t-1,ξ,υ。

由于一定程度上从时间序列间接观测隐含变量，由于x_t|t-1、测量函数g和测量噪声的不确定性，第a-维输出仍然是随机的。然而，只有不确定输入是高斯分布的，预测测量可以解析地计算[23]。假设p(x_t|t-1|y_1:t-1)和p(y_t|x_t|t-1,ξ,υ)是高斯分布的，那么p(y_t|y_1:t-1,x_t|t-1,ξ,υ)＝∫p(y_t|x_t|t-1,ξ,υ)p(x_t|t-1|y_1:t-1)dx_t|t-1仍然是高斯分布。因此p(y_t|y_1:t-1,x_t|t-1,ξ,υ)的预测均值可以计算如下[31]

${\mathrm{\μ}}_{t|t-1}^{y_a}{=E}_{x_{t|t-1,g}}\left(g\left(x_{t|t-1}\right)\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x,\mathrm{\ξ},\mathrm{\υ})={\left({(K_g^{\left(a\right)}+{\mathrm{\σ}}_v^{2}I)}^{-1}{\mathrm{\υ}}^{\left(a\right)}\right)}^T{q}_g^{\left(a\right)}---\left(24\right)$

这里υ^(a)是υ(a＝1,2,...,E)的第a-列，ξ_i是ξ的第i-行。 $q_g^{\left(a\right)}=[q_1^{\left(a\right)g},q_2^{\left(a\right)g},L,q_N^{\left(a\right)g}],$ 其i^th元素为

$q_i^{\left(a\right)g}=E_{x_{t|t-1}}\left(k_{{\mathrm{\α}}_a}(x_{t|t-1},{\mathrm{\ξ}}_i)\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)={\mathrm{\α}}_a^2{|{\mathrm{\Σ}}_{t|t-1}^x{\mathrm{\Λ}}_a+I|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{({\mathrm{\ξ}}_i-{\mathrm{\μ}}_{t|t-1}^x)}^T{({\mathrm{\Σ}}_{t|t-1}^x+{\mathrm{\Λ}}_a)}^{-1}({\mathrm{\ξ}}_i-{\mathrm{\μ}}_{t|t-1}^x)).$

注意到测量函数与测量噪声相互独立的，与之间的协方差可以通过如下计算得到

$\begin{array}{c}{\mathrm{\Σ}}_{t|t-1}^{y_a{y}_b}=\mathrm{cov}(y_{t|t-1}^{\left(a\right)},y_{t|t-1}^{\left(b\right)}|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)=\mathrm{cov}(g_{\left(a\right)}\left(x_{t|t-1}\right)+w_a,g_{\left(b\right)}\left(x_{t|t-1}\right)+w_b|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)\\ =\mathrm{cov}(g_{\left(a\right)}\left(x_{t|t-1}\right),g_{\left(b\right)}\left(x_{t|t-1}\right)|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)-{\mathrm{\μ}}_{t|t-1}^{y\left(a\right)}{\mathrm{\μ}}_{t|t-1}^{y\left(b\right)}+{\mathrm{\δ}}_{a,b}{\mathrm{\σ}}_w^{2}I\end{array},---\left(25\right)$

y_t|t-1对应的协方差矩阵为

${\mathrm{\Σ}}_{t|t-1}^y=\left[\begin{array}{ccc}{\mathrm{\Σ}}_{t|t-1}^{y_1}& L& {\mathrm{\Σ}}_{t|t-1}^{y_1{y}_E}\\ M& O& M\\ {\mathrm{\Σ}}_{t|t-1}^{y_E{y}_1}& L& {\mathrm{\Σ}}_{t|t-1}^{y_E}\end{array}\right]$

为简洁起见，g_(a)(x_t|t-1)记为g_(a)(x_t|t-1)和g_(b)(x_t|t-1)协方差为

$\begin{array}{c}\mathrm{cov}(g_{t|t-1}^{\left(a\right)},g_{t|t-1}^{\left(b\right)}|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)=E_{x_{t|t-1},g}\left(g_{t|t-1}^{\left(a\right)}{g}_{t|t-1}^{\left(b\right)}\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)\\ =E_{x_{t|t-1}}\left\{E_g\left(g_{t|t-1}^{\left(a\right)}\right|x_{t|t-1})E_g\left(g_{t|t-1}^{\left(b\right)}\right|x_{t|t-1})\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x\}\end{array}---\left(26\right)$

考虑到

$E_g\left(g_{t|t-1}^{\left(a\right)}\right|x_{t|t-1})=\∫p\left(g\right)\left(g_{t|t-1}^{\left(a\right)}\right|x_{t|t-1})\mathrm{dg}=[k_{{\mathrm{\α}}_a}(x_{t|t-1},{\mathrm{\ξ}}_i),L,k_{{\mathrm{\α}}_a}(x_{t|t-1},{\mathrm{\ξ}}_M)]{(K_g^{\left(a\right)}+{\mathrm{\σ}}_v^{2}I)}^{-1}{\mathrm{\υ}}^{\left(a\right)}---\left(27\right)$

令 $Q_{\mathrm{ij}}^{g(a,b)}={\left[q_{i,j}^{g(a,b)}\right]}_{i,j=\mathrm{1,2}}^N,{\hat{z}}_{\mathrm{ij}}={\mathrm{\Λ}}_b{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}{\mathrm{\ξ}}_i+{\mathrm{\Λ}}_a{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}{\mathrm{\ξ}}_j,$ 我们有

$\begin{array}{c}{q}_{i,j}^{g(a,b)}=E_{x_{t|t-1}}\left(k_{{\mathrm{\α}}_a}(x_{t|t-1},{\mathrm{\ξ}}_i)k_{{\mathrm{\α}}_b}(x_{t|t-1},{\mathrm{\ξ}}_j)\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)\\ ={\mathrm{\α}}_a^2{\mathrm{\α}}_b^2{|({\mathrm{\Λ}}_a^{-1}+{\mathrm{\Λ}}_b^{-1}){\mathrm{\Σ}}_{t|t-1}^x+I|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{({\mathrm{\ξ}}_i-{\mathrm{\ξ}}_j)}^T{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}({\mathrm{\ξ}}_i-{\mathrm{\ξ}}_j))\\ \×\exp (-\frac{1}{2}{({\hat{z}}_{\mathrm{ij}}-{\mathrm{\μ}}_{t|t-1}^x)}^T{({({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}+{\mathrm{\Σ}}_{t|t-1}^x)}^{-1}({\hat{z}}_{\mathrm{ij}}-{\mathrm{\μ}}_{t|t-1}^x))\end{array}---\left(28\right)$

因此，和的协方差写为

$\mathrm{cov}(g_{t|t-1}^{\left(a\right)},g_{t|t-1}^{\left(b\right)}|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)={\left({(K_g^{\left(a\right)}+{\mathrm{\σ}}_g^{2}I)}^{-1}{\mathrm{\υ}}^{\left(a\right)}\right)}^T{q}_{i,j}^{g(a,b)}\left({(K_g^{\left(b\right)}+{\mathrm{\σ}}_g^{2}I)}^{-1}{\mathrm{\υ}}^{\left(b\right)}\right)---\left(29\right)$

根据上面结论，p(y_t|y_1:t-1,x_t|t-1,ξ,υ)的和就可以精确得到。

3.2状态滤波

滤波的目的就是使用当前和过去的观测逼近后验分布p(x_t|y_1:t)以减少隐藏状态x_t的不确定性。后验分布p(x_t|y_1:t)可以由Bayesian定理计算

$p\left(x_t\right|y_{1:t})=\frac{p\left(y_t\right|x_t)p\left(x_t\right|y_{1:t-1})}{p\left(y_t\right|y_{1:t-1})}=\frac{p(y_t,x_t|y_{1:t-1})}{p\left(y_t\right|y_{1:t-1})}---\left(30\right)$

这里p(y_t|x_t)由式Eq。(21)计算得到，p(x_t|y_1:t-1)在状态预测步得。注意到p(y_t|y_1:t-1)＝∫p(x_t|y_1:t-1)p(y_t|x_t)dx_t是归一化项并且(22)的分子不能精确计算，式(22)不能得到闭解，因此使用高斯分布逼近p(x_t|y_1:t)。

p(x_t|y_1:t)的参数由下式计算

${\mathrm{\μ}}_{t|t}^x={\mathrm{\μ}}_{t|t-1}^x+{\mathrm{\Σ}}_{t|t-1}^{\mathrm{xy}}{\left({\mathrm{\Σ}}_{t|t-1}^y\right)}^{-1}(y_t-{\mathrm{\μ}}_{t|t-1}^y)---\left(31\right)$

${\mathrm{\Σ}}_{t|t}^x={\mathrm{\Σ}}_{t|t-1}^x+{\mathrm{\Σ}}_{t|t-1}^{\mathrm{xy}}{\left({\mathrm{\Σ}}_{t|t-1}^y\right)}^{-1}{\mathrm{\Σ}}_{t|t-1}^{\mathrm{yx}}---\left(32\right)$

滤波算法计算状态的均值、协方差和互协方差(cross-covariance)以及测量预测分布。和在测量预测中得到。下面计算状态预测分布的互协方差均值和协方差

从y_1:t-1可以计算状态和测量预测在t时刻的互协方差

${\mathrm{\Σ}}_{t|t-1}^{\mathrm{xy}}=[\mathrm{cov}(x_{t|t-1},g_1\left(x_{t|t-1}\right)+w_1),L,\mathrm{cov}(x_{t|t-1},g_E\left(x_{t|t-1}\right)+w_E)|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x],$ 即

$\mathrm{cov}(x_{t|t-1},g_a\left(x_{t|t-1}\right)+w_a|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)=E_{x_{t|t-1},g}(x_{t|t-1},g_a\left(x_{t|t-1}\right)|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x)-{\mathrm{\μ}}_{t|t-1}^x{\mathrm{\μ}}_{t|t-1}^{y\left(a\right)}$

$E_{x_{t|t-1},g_a}\left[x_{t|t-1}{g}_a\left(x_{t|t-1}\right)\right|{\mathrm{\μ}}_{t|t-1}^x,{\mathrm{\Σ}}_{t|t-1}^x]=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\τ}}_i^{\left(a\right)g}{q}_i^{\left(a\right)g}{\mathrm{\Σ}}_{t|t-1}^x{({\mathrm{\Σ}}_{t|t-1}^x+{\mathrm{\Λ}}_a)}^{-1}({\mathrm{\ξ}}_i-{\mathrm{\μ}}_{t|t-1}^x)\∈R^D---\left(33\right)$

这里 ${\mathrm{\τ}}_a={\left[{\mathrm{\τ}}_i^{\left(a\right)g}\right]}_{i=\mathrm{1,2}}^N={\left(K_g^{\left(a\right)}\right)}^{-1}{\mathrm{\υ}}^{\left(a\right)}.$

假定FITCSGP回归预测辨识系统模型，记转换模型的FTICGP回归超超参数为θ_f，相应的伪输入训练数据为α,β∈R^N×D以及FITCGP为p(x_t+1|y_1:t,α,β)。在时间t出的状态预测为

$p\left(x_t\right|y_{1:t-1},\mathrm{\α},\mathrm{\β})=p\left(x_{t-1}\right|y_{1:t-1})p\left(x_t\right|x_{t-1},\mathrm{\α},\mathrm{\β}){\mathrm{dx}}_{t-1}=E_{x_{t-1}}\left(p\left(x_t\right|x_{t-1},\mathrm{\α},\mathrm{\β})\right)---\left(34\right)$

显然，p(x_t|x_t-1)是由转换函数计算的转换概率，p(x_t-1|y_1:t-1)是逼近高斯函数，他们的积分不是高斯分布。这里让人用高斯分布逼近该概率

对SGP，伪训练数据能够被解释为N对来自x_t-1→x_t的转换函数的独立的观测对，以及N对来自x_t→y_t的测量函数的独立观测。并且带有伪训练集的稀疏GPS可以施加条件独立x_t+1⊥x_t-1|x_t,α,β和a,b＝1,2,L,D。这里f_a(x_t)是f(x_t)的第a-为元素。给出考虑到转换函数和状态x_t-1的不确定性，f(x_t-1)的每一维均值，单独计算为

${\mathrm{\μ}}_{t|t-1}^{x\left(a\right)}=E_{x_{t-1,f}}\left(f_a\left(x_{t-1}\right)\right|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})---\left(35\right)$

相似地，x_t的协方差为

${\mathrm{\Σ}}_{t|t-1}^x\left|\right|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1,t-1}^x,\mathrm{\α},\mathrm{\β}=\left[\begin{array}{ccc}\mathrm{var}\left(x_t^{\left(1\right)}\right|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})& L& \mathrm{cov}(x_t^{\left(1\right)},x_t^{\left(D\right)}{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})\\ & O& M\\ \mathrm{cov}(x_t^{\left(1\right)},x_t^{\left(D\right)}|{\mathrm{\μ}}_{t-1|t-1}^{\left(D\right)},{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})& L& \mathrm{var}\left(x_t^{\left(D\right)}\right|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})\end{array}\right]$

$\begin{array}{c}\mathrm{cov}(x_t^{\left(x\right)},x_t^{\left(b\right)}|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})=\mathrm{cov}(f_a\left(x_{t-1}\right)+w_a,f_b\left(x_{t-1}\right)+w_b|{\mathrm{\μ}}_{t|t}^x,{\mathrm{\Σ}}_{t|t}^x,\mathrm{\α},\mathrm{\β})\\ =\mathrm{cov}(f_a\left(x_{t-1}\right),f_b\left(x_{t-1}\right)|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x,\mathrm{\α},\mathrm{\β})-{\mathrm{\μ}}_{t|t-1}^{x\left(a\right)}{\mathrm{\μ}}_{t|t-1}^{x\left(b\right)}+{\mathrm{\δ}}_{a,b}{\mathrm{\σ}}_w^2\end{array}---\left(36\right)$

和形式类似3。1部分，为简单起见，其详细公式省略。记为在时间t的状态预测协方差由如下方式计算：

$\begin{array}{c}\mathrm{cov}(f_{t|t-1}^{\left(a\right)},f_{t|t-1}^{\left(b\right)}|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x)=E_{x_{t-1|t-1},f}(f_{t|t-1}^{\left(a\right)}{f}_{t|t-1}^{\left(b\right)}\text{|}{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x)\\ =E_{x_{t-1|t-1}}\left\{E_f\left(f_{t|t-1}^{\left(a\right)}\right|x_{t-1|t-1})E_f\left(f_{t|t-1}^{\left(b\right)}\right|x_{t-1|t-1})\right|{\mathrm{\μ}}_{t-1|t-1}^x,{\mathrm{\Σ}}_{t-1|t-1}^x\}\\ ={\left({(K_f^{\left(a\right)}+{\mathrm{\σ}}_w^{2}I)}^{-1}{\mathrm{\β}}^{\left(a\right)}\right)}^T{q}_{i,j}^{(a,b)f}\left({(K_f^{\left(b\right)}+{\mathrm{\σ}}_w^{2}I)}^{-1}{\mathrm{\β}}^{\left(b\right)}\right)\end{array}---\left(37\right)$

其中

$E_f\left(f_{t|t-1}^{\left(a\right)}\right|x_{t-1|t-1})=[k_{{\mathrm{\β}}_a}(x_{t-1|t-1},a_1),L,k_{{\mathrm{\β}}_a}(x_{t-1|t-1},a_N)]{(K_f^{\left(a\right)}+{\mathrm{\σ}}_w^{2}I)}^{-1}{\mathrm{\β}}^{\left(a\right)}$

$\begin{array}{c}{q}_{i,j}^{(a,b)f}=E_{x_{t-1|t-1}}\left(k_{{\mathrm{\α}}_a}(x_{t-1|t-1},a_i)k_{{\mathrm{\α}}_b}(x_{t-1|t-1},a_j)\right|{\mathrm{\μ}}_{t|t}^x,{\mathrm{\Σ}}_{t|t}^x)\\ ={\mathrm{\α}}_a^2{\mathrm{\α}}_b^2{|({\mathrm{\Λ}}_a^{-1}+{\mathrm{\Λ}}_b^{-1}){\mathrm{\Σ}}_{t-1|t-1}^x+I|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{(a_i-a_j)}^T{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}(a_i-a_j))\\ \×\exp (-\frac{1}{2}{({\hat{z}}_{\mathrm{ij}}-{\mathrm{\μ}}_{t-1|t-1}^x)}^T{({({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}+{\mathrm{\Σ}}_{t-1|t-1}^x)}^{-1}({\hat{z}}_{\mathrm{ij}}-{\mathrm{\μ}}_{t-1|t-1}^x))\end{array}$

${\hat{z}}_{\mathrm{ij}}={\mathrm{\Λ}}_b{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}{\mathrm{\α}}_i+{\mathrm{\Λ}}_a{({\mathrm{\Λ}}_a+{\mathrm{\Λ}}_b)}^{-1}{\mathrm{\α}}_j.$

基于上述结果，滤波后状态x_t|t可以从(31)和(32)给出当前观测y_t。

2.3实现细节

SGP学习的初始隐含状态：

算法使用2个FITCSGP回归模型对状态空间的系统和测量模型进行建模并同时优化SGP超参数和伪输入。状态并不能直接从时间序列观测到，因此有可能陷入局部最优解。当初始值在最优解的域内才能得到好的解。因此合理的隐藏状态对学习GPIL非常关键。这里使用线性状态空间辨识隐含状态和参数作为GPIL初始值是一个较好的选择。N4SID算法被认为一种流形的子空间辨识方法(SIM)，该方法基于SVD方法辨识状态空间模型的闭解形式。N4SID算法能充分接近真实隐藏状态，在很多实际应用作为非线性动态系统的初始方法。线性状态空间训练后，使用Kalman滤波方法估计隐藏状态。

SGPs超参数初始化通常使用训练数据样本子集使用GP拟合得到，这样就可以减少计算代价。由于隐藏状态可以看做不确定输入，标准SGP学习只是优化伪输入，降低了过拟合可能性。然后通过估计SGP在伪输入出的伪输出均值。一旦使用2个SGP回归模型估计出系统模型和测量模型，可以交替置信测量和滤波步骤使用GP-ADF算法从测量数据中估计隐藏状态。

GPIL算法的学习过程总结如下:

输入：测量时间序列Y∈R^T×E，隐藏状态维数D，伪样本数量M，终止阈值ε

输出:延长状态X∈R^T×E，对应系统和测量函数f(·)andg(·)的SGPs回归模型超参数θ_f和θ_g伪输入和输出分别为{ξ,ν}及{α,β}。

1)初始化隐藏状态

1.1)使用N4SID算法计算M维初始隐藏状态和参数；

1.2)使用Kalman滤波对初始状态进行滤波。

2)训练2个SP回归模型得到GPADF的初始值

2.1)根据隐藏状态x_t学习SGP估计系统模型，所学超参数和伪输入、输出分别记为θ_f和{ξ,ν}；

2.2)根据隐藏状态x_t和测量y_t学习SGP估计测量模型，学习的超参数和伪输入输出分别记为θ_g和{α,β}。

3)使用梯度下降方法，利用上面的初始值对最大化边缘似然函数 $L\left(\mathrm{\Θ}\right)=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\log p\left(y_t\right|y_{1:t-1},\mathrm{\Θ})=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\log N\left(y_t\right|{\mathrm{\μ}}_{t|t-1}^y,{\mathrm{\Σ}}_{t|t-1}^y)$ 进行优化求解最优超参数，n＝1

重复

4)运行GP-ADF进行状态滤波

4.1)从事先给定均值和协方差的高斯分布产生初始状态x₁；

4.2)t＝1到t＝T；

a)使用参数和伪输入、输出的θ_g,{α,β}SGP，使用式(24)，(25)and(29)，(33)所示根据x_t预测测量y_t，计算时间t时刻预测测控的均值和协方差

b)使用(31)和(32)从y_1:t对状态x_t滤波；

c)使用参数θ_f和伪输入、输出{ξ,ν}SGP预测状态从x_t到x_t+1，使用(35)和(36)计算t+1时刻的状态预测均值和协方差

${\mathrm{\Θ}}^{(n+1)}={\mathrm{\Θ}}^{\left(n\right)}-\mathrm{\λ}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\left\{\right[\frac{1}{2}{(y_t-{\mathrm{\μ}}_{t|t-1}^y)}^T{\left({\mathrm{\Σ}}_{t|t-1}^y\right)}^{-1}\frac{{\∂\mathrm{\Σ}}_{t|t-1}^y}{{\∂\mathrm{\Θ}}^{\left(n\right)}}+\frac{{\∂\mathrm{\μ}}_{t|t-1}^y}{{\∂\mathrm{\Θ}}^{\left(n\right)}}]{\left({\mathrm{\Σ}}_{t|t-1}^y\right)}^{-1}(y_t-{\mathrm{\μ}}_{t|t-1}^y)-\frac{1}{2}\mathrm{tr}\left({\left({\mathrm{\Σ}}_{t|t-1}^y\right)}^{-1}\frac{{\∂\mathrm{\Σ}}_{t|t-1}^y}{{\∂\mathrm{\Θ}}^{\left(n\right)}}\right)$

n＝n+1

直到||Θ^(n-1)-Θ⁽ⁿ⁾||≤ε

5)返回X，SGPs超参数和伪训练数据Θ。

一旦最优SGPs超参数和伪训练数据获取，使用测量预测和状态滤波步骤新的状态很容易从当前输出测量y_T+1计算x_T+1。

Claims (1)

1.一种基于矿用锂离子电池的电源管理方法，其特征是：所述的电源管理方法：主控单元和单总线温度采集单元，在电源管理系统首次上电后，主控单元加电自检并配置电源管理单元的参数，参数配置完成后，系统等待具有唯一序列号的单总线温度传感器接入，主控单元间隔搜索单总线上的温度传感器数量及各自的序列号，并存储在具有掉电不丢失数据特性的FLASH中，此处是由主控单元的FLASH模拟EEPROM实现；

电流采集单元，根据采集的数据和系统以及不同电流等级和精度要求，共分为两组不同的电流采集方案：a、充电电流采集采用霍尔电流采集隔离方案；b、放电和容量自检电流采样通过串联精密电阻和电压变换测量，电流采集是基于主控单元自带的12位AD模块，外部加载独立高精度的电压源，结合特定的采样方案和数据处理方法，实现高精度的电流采样；

所述的电源管理单元，主控单元通过SPI接口操作电源管理单元，实现电压采集和充电均衡管理任务，经任务调度，可以不停的回读单体电池的电压、电池组的总电压；充电时，根据单体电池不同的电压等级，主控单元选择不同的均衡阀值来进行充电均衡管理，同时主控单元判断各部分的电压采集是否正常，以此判断电压采集线是否出问题，并更新温度和状态数据；

所述的数据传输单元，通过RS232接口与主控单元进行通信，接收电池数据，用于数据的显示，并给主控单元发出参数设置指令，可以设置调整的参数有：单体电池过压阀值、单体电池欠压阀值、充放电电流阀值、充电均衡时间、过流保护恢复时间、充电均衡电压、电池组容量值、充放电和容量自检电流修正系数。