线性回归方法RC电路的阶跃响应—时间常数τ

突然有个想法用机器学习的基本方法——线性回归法来学习一阶RC电路的阶跃响应，从而得到RC电路的结构特征——时间常数τ(即 R*C)。答案无疑是肯定的，但问题是如何通过最小二乘法和正规方程来减少信号采集噪声对更多采样点的τ估计值的影响。另外，由于最近对 Jupyter 和 numpy 的摆弄，我只是尝试使用 Jupyter 而不是 matlab。结果发现它非常容易使用，尤其是能够在 Jupyter 脚本中插入 LaTeX 格式的公式，非常方便！我尝试了一个常用的工具软件，直接将纸上的手写公式转换成 Jupyter 脚本。

以下原创内容欢迎网友转载，但请注明出处：

一、理论推导

1。线性回归分析和正态方程

传统上，线性回归问题是使用最小二乘法（即正规方程组）来解决最小化线性方程组的均方误差的问题。已知输出输入X为由多个变量组成的行向量，W为系数向量（列向量），b为偏置

在机器学习中，将每个输入 x 作为一行，形成一个更大的矩阵偏最小二乘回归的线性与非线性方法，其中每一行代表一个样本，称为设计矩阵 X(train)。如果样本数为k，那么X(train)有k行，每个样本会得到一个输出y，y被收集成一个列向量Y(train)，k个相同的bs也组成一个列向量b。为了简化表达式，b被简化为一个常数b附加到系数列向量W的最后一个，X(train)在每行的末尾添加一个1，W(column vector)在每个行的末尾添加一个待估计的b结尾。为了最小化估计结果之间差异的平方和：

和Y(train)，有正规方程可以求解W：

2.一阶RC电路跳转响应的顺序

一阶RC电路的电路模型如下图所示。

图。 1 一阶RC电路

幅度为Vcc的阶跃信号从Vin输入，输出在Vout处测量。求解微分方程得到一个响应函数，其参数是时间 t。

其中时间常数 τ = R*C。我想通过在阶跃信号输入条件下测量实际 RC 电路的响应曲线 V(t) 并通过 V(t) 来估计时间常数 τ。如果做纯理论推导，只要知道t0时刻的输出电压V(t0)，就可以解方程(2）得到τ。但在实际电路中，电压 V(t0) ) 测量往往受到很多干扰，不准确，是否可以测量多个 t 值矩对应的 V(t)偏最小二乘回归的线性与非线性方法，使用正规方程 (1）获得统计上的最优估计

怎么样？这是下一个要解决的问题。

3.非线性函数的最小二乘估计

仔细观察目标函数（0）公式适用于正规方程的特性，可以发现如果要制作非线性公式（2）公式可以使用正规方程，你必须让：

1）包含待估计变量τ的函数作为(0）公式中的系数W，而设计矩阵X(train)可以作为包含时间t的函数或测量电压V(t)。

2）W和X(train)的关系一定是简单的乘法。

显然，只有时间 t 的序列作为设计矩阵 X(train)，有可能 W 和 X(train) 之间的关系必须相乘。至于其他非线性部分，可以通过等价变换变换到方程的另一边，作为Y(train)。综上所述，(2）公式可以转化为(3）公式。

(整个3）公式左边可以计算得到Y(train)；(<@公式3）右边时间t的序列形成设计矩阵X( train)，1/τ形成系数矩阵W等价于(0）公式。这样，正规方程(2）公式求解统计意义上的最优估计

。当然，从(3）公式中也可以看出，经过线性校正后，系数向量W在模型中只有一个元素，是一个标量，偏置b也总是等于0。

二、仿真模型

想用最近尝试用Jupyter和numpy替换熟悉的Matlab来验证上面非线性函数最小二乘估计的思想。让我们先建立一个模型：

1） 输入是Vcc的大小是3.3V

2）时间常数τ为0.1秒；

3）简单模拟的采样间隔随机性：是的间隔等于delta1=0.0015秒和delta1=0.0011秒这两个序列的叠加。

4）采样总长度为n=5乘以τ；

5） 信号中叠加了幅度为20mV左右的白噪声——至于为什么是20mV，后面会详细讲解。

使用python和numpy进行数值模拟的代码如下：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 tao=0.1#时间常数
 4 vcc=3.3#电源电压
 5 n=5#时长取时间常数tao的n倍
 6 delta1=0.0015#第一种采样间隔
 7 delta2=0.0011#第一种采样间隔
 8 t1=delta1*np.arange(np.ceil(n*tao/delta1))
 9 t2=delta2*np.arange(np.ceil(n*tao/delta2))
10 t=np.append(t1,t2)#为演示最小二乘拟合的功能，故意结合两种采样率下的时间点
11 t.sort()#对t进行排序
12 plt.plot(t)

13 s=vcc*(1-np.exp(-t/tao))#理论的波形曲线
14 plt.plot(t,s)#注意这里的plot函数使用了x轴和y轴两个轴，因为s中的数据不是均匀的
15 N_amp=np.exp(-n)*vcc
16 N_amp
17 noise = np.random.uniform(-N_amp, N_amp, (len(t)))#噪声：正负np.exp(-5)*3.3之间均匀分布
18 s_nr=s+noise#加入噪声后的信号
19 plt.plot(t,s_nr)
20 yt=np.log(vcc/(vcc-s_nr))
21 plt.plot(t,yt)
22 yt=np.mat(yt)
23 yt=yt.T
24 x=np.mat(t)#将X转换为矩阵数据类型
25 x=x.T#正规方程中的x应该是个列向量
26 w=(np.linalg.inv(x.T*x))*x.T*yt#求解正规方程
27 E_tao = np.round(1/float(w),4)#对时间常数的tao的估计，保留到4位小数
28 E_tao

非线性函数的最小二乘拟合

说明：

1） 序列包含两个等差数列t1和t2的融合，它们的区间互质，不重复，意在模拟采样时间的随机性。最后，使用 sort() 方法对时间序列进行排序的目的是为了更容易更直观地观察波形。如果只是使用正规方程，则不需要重新排序。

2） 修正后的非线性方程（3）与原方程（2））相比有一个重大缺陷：值不能为负——如果是纯理论推导，这不能但是如果随机噪声后的V(t)有可能大于步幅Vcc，那么括号内就会变成负数，这样计算就变得毫无意义了。我的解决办法是限制假设随机的幅度噪声对仿真截止时间V(t)和Vcc的差值，以及代码中的N_amp。由于仿真结束时间是n(= 5）tau，所以N_amp的值等于np.exp( -n)*vcc。

这没有理论依据，只受限于(3）和(2）方程之间的非完全等价变换——线性化校正的代价。

3）由于只有一个参数需要估计（1/τ），所以正规方程的解也是一个只有一个元素的矩阵。将其转换为标量并将其反转以获得最佳估计

.

三、结果评估

加噪声后的信号如下图所示，与通常情况下的实测波形吻合较好。

图 2 仿真产生的噪声 的阶跃响应

对这些添加的噪声信号进行线性校正后，通过最小二乘估计得到的信号yt如下图所示。正如我所料，基本趋势是一条斜率为 (1/τ) 的直线。

图3 图2线性校正后待估计的信号

但是可以清楚的看出，由于(3）公式左边在V(t)的大小接近Vcc时，会放大加入的白噪声。因此，当t增大时，由白噪声引起的信号大大增加。也就是说，当应用正规方程时，t 值较大时噪声的影响要大于 t 值较小时的影响。这是使用 a 付出的重要代价非线性校正函数。

另外，多次运行上述代码，得到的值都比实际的τ略小（=0.1）—大约0.099左右，即, 不是无偏估计。这应该是由于线性校正函数((3）公式左边),是对噪声大于0小于0的部分进行了不对称变换造成的. 虽然由此造成的偏差不大，但这也是使用非线性校正功能需要付出的代价。

四、Jupyter 笔记本

上述练习的一个重要目的是练习使用 Jupyter 笔记本并嵌入交互性良好的公式等信息。下面是部分程序的运行效果截图。虽然对markdown语法不熟悉，格式和布局不够美观，但明显可以提高代码的可读性。

图4 Jupyter notebook运行效果截图

需要记录的是公式代码的嵌入过程：

1）我第一次写论文我手写了一些公式，用手机拍了照片，比如：

图 5 手写公式

2）使用 mathpix 工具对这些照片进行截图并扫描（mathpix 工具适用于 Windows 和 iOS，两者都可以免费试用）。 Mathpix 可以直接获取 LaTeX 格式的公式表达式。下图是iOS版mathpix界面截图。

图6 iOS版mathpix截图

3）在Jupyter notebook顶部的下拉菜单中选择单元格格式为Markdown；将LaTeX格式的公式表达式粘贴到单元格中，在LaTeX公式表达式前后添加“$$”符号，运行单元格即可得到公式，效果如图4所示。

图 7 在 Jupyter notebook 中输入 LaTeX 公式

五、进一步的实际测试

工作还没有结束。应该测试一个简单的实际电路。后续假期抽半天在STM32开发板上完成这个工作：用GPIO生成一个保存信号后，连续采集5个τ时间长度的信号，代入正规方程求解，欢迎大家继续关注更新吧。

…