锂电池研究之七——基于 Pytorch 的高斯函数拟合时间序列数据

原创博客，转载请注明出处，谢谢！

代码下载地址：https://github.com/XiuzeZhou/CALCE

1. 需求

最近有同学发邮件问我：因为我不是计算机方向的，不熟悉机器学习的方法，所以想用传统的数据拟合的方法来对锂电池的数据进行建模——输入是时间变量，输出是电池容量，用 python 该怎么去实现？

这个问题，不是很简单吗？传统的数据拟合应该很多人研究才对，Pytorch 这么强大，做数据拟合不是大炮打蚊子吗，网上找找应该到处都是呀。然后我发现，过于乐观了。上网一搜，发现网上的数据拟合大部分是 matlab 的代码，少部分有 python 的代码不是跑不通，就是和时间序列数据没半毛钱关系，更不要说 Pytorch 代码了！

于是，抽空花半天时间写了一个用 Pytorch 来计算高斯函数，拟合时间序列数据的项目。

2. 问题定义

根据该同学提供的论文，其中的高斯函数定义如下：

该高斯函数由 3 组正态分布构成（仅用正态分布表示高斯函数 f(x) 中的每一项高斯分布，与定义区别开而已）。其中，$a_i,b_i,c_i,(i=1,2,3)$ 是三组需要拟合的参数，x 则是时间变量，在锂电池中又可以叫充放电次数或者等效使用时间。

3. 时间序列数据

时间序列数据，是随着时间变化的一系列数据。在锂电池中，是其随着充放电次数的增加性能衰减的过程。本篇博文，以 CALCE 锂电池数据作为说明。CALCE 原始数据集的数据格式是 xlsx 的，数据处理部分就不在赘述了，想了解这这部分内容的同学，请看我的博客： https://snailwish.com/395/。

最后，我们可以得到锂电池的电池容量随着充放电次数的衰减趋势图如下所示：

其中，横坐标是充放电次数（t），纵坐标是电池容量。

4. 代码及结果

I. 网络定义

MLP 的定义很简单，主要就是构建前面的高斯拟合函数。

class Net(nn.Module):
    def __init__(self, k=3, eps=1e-6):
        super(Net, self).__init__()
        self.eps = eps
        self.w = torch.nn.Parameter(torch.randn((3, k)), requires_grad=True)

    def forward(self, x): 
        out = 0
        for i in range(k):
            out += self.w[0, i] * torch.exp(
                -torch.pow((x - self.w[1, i])/(self.w[2, i] + self.eps), 2)) 
        return out

其中，k 表示高斯函数中有 3 个正态分布函数，eps 表示一个极小的数字，用来避免分母为零的情况，w 表示需要拟合的参数, forward(x) 中定义的是前面的高斯函数 f(x)。

II. 训练模型

def train(data, k=3, lr=1e-4, stop=1e-3, epochs=100000, device=device):  
    model = Net(k=k)
    model = model.to(device)
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    criterion = nn.MSELoss()

    loss_list = []
    epoch = 0
    while True:
        train_data = np.random.permutation(data)
        x, y = train_data[:,0], train_data[:,1]
        X = np.reshape(x, (-1, 1)).astype(np.float32)
        y = np.reshape(y, (-1, 1)).astype(np.float32)
        X = Variable(torch.from_numpy(X)).to(device)
        y = Variable(torch.from_numpy(y)).to(device)
        y_ = model(X)
        loss = criterion(y_, y)
        #print(loss.detach().numpy())
        optimizer.zero_grad()          # clear gradients for this training step
        loss.backward()                # backpropagation, compute gradients
        optimizer.step()               # apply gradients

        if epoch % 100 == 0:
            loss_list.append(loss.detach().numpy())

        if (loss.detach().numpy() < stop) or (epoch > epochs):
            break

        epoch +=1

    return model, loss_list

其中，k 表示高斯函数中有 3 个正态分布函数，lr 表示学习率，stop 表示损失的停止条件，epochs 表示运行的最大次数，device 表示是设备 CPU 还是 GPU。

model, loss_list = train(
    data=np.c_[x, y], k=k, lr=lr, stop=stop, epochs=epochs, device=device)
print('Optimal parameters: ', model.w)
>>
a: [ 106.4594,  -91.1130,  -19.7342]
b: [ 286.9670,  314.5258,  152.2022]
c: [ 302.1231,  289.4315, -262.6197]

model.w 的结果就是求得的 a, b 和 c 的参数值。注意，你的结果并不一定就是一模一样的值，运行多次，可能会有不同的值。准确点数学上描述就是求值的极值点可能存在多个。

III. 预测效果

def predict(model, x, device=device):
    x = np.array(x).astype(np.float32)
    x= torch.from_numpy(x).to(device)
    model.eval()
    torch.no_grad()
    y_ = model(x)
    return y_.detach().numpy()

我们先看一下训练过程中的损失函数的变化，如下图所示：

由曲线可知，损失函数确实收敛了，说明这些参数训练和求值没什么问题。然后，看一下最后用这些参数拟合的效果：

预测曲线光滑，基本趋势拟合的很不错！