突然有个想法, 利用机器学习的基本方法 -- 线性回归方法, 来学习一阶 RC 电路的阶跃响应, 从而得到 RC 电路的结构特征 -- 时间常数τ(即 R*C). 回答无疑是肯定的, 但问题是怎样通过最小二乘法, 正规方程, 以更多的采样点数来降低信号采集噪声对τ估计值的影响. 另外, 由于最近在捣鼓 Jupyter 和 numpy 这些东西, 正好尝试不用 matlab 而用 Jupyter 试试看. 结果是意外的好用, 尤其是在 Jupyter 脚本中插入 LaTeX 格式的公式的功能, 真是太方便了! 尝试了直接把纸上手写的公式转换到 Jupyter 脚本中的常见工具软件.
一, 理论推导
1.线性回归分析及正规方程
传统意义说, 线性回归问题是用最小二乘法(即正规方程), 解决线性方程组的均方误差最小化问题. 已知输出输入 X 是由多个变量构成的行向量, W 是系数向量(列向量),b 为偏置
在机器学习中, 把每次的输入 x 作为一行组成更大的矩阵, 即每一行代表一个样本, 该矩阵称为设计矩阵 X(train). 若样本数为 k, 则 X(train)有 k 行, 每个样本都会得到一个输出 y, 将 y 集合成一个列向量 Y(train),k 个相同的 b 也组成列向量 b. 为简化表达, 将 b 简化为附加在系数列向量 W 最后的常数 b,X(train)则在每行的最后增加一个 1,W(列向量)的最后增加一个待估计的 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/τ则构成相当于 (0) 式中的系数矩阵 W. 这样就可以通过正规方程 (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 进行数值仿真的代码如下:
- import numpy as np
- import matplotlib.pyplot as plt
- tao=0.1# 时间常数
- vcc=3.3# 电源电压
- n=5# 时长取时间常数 tao 的 n 倍
- delta1=0.0015# 第一种采样间隔
- delta2=0.0011# 第一种采样间隔
- t1=delta1*np.arange(np.ceil(n*tao/delta1))
- t2=delta2*np.arange(np.ceil(n*tao/delta2))
- t=np.append(t1,t2)# 为演示最小二乘拟合的功能, 故意结合两种采样率下的时间点
- t.sort()# 对 t 进行排序
- plt.plot(t)
- s=vcc*(1-np.exp(-t/tao))# 理论的波形曲线
- plt.plot(t,s)# 注意这里的 plot 函数使用了 x 轴和 y 轴两个轴, 因为 s 中的数据不是均匀的
- N_amp=np.exp(-n)*vcc
- N_amp
- noise = np.random.uniform(-N_amp, N_amp, (len(t)))# 噪声: 正负 np.exp(-5)*3.3 之间均匀分布
- s_nr=s+noise# 加入噪声后的信号
- plt.plot(t,s_nr)
- yt=np.log(vcc/(vcc-s_nr))
- plt.plot(t,yt)
- yt=np.mat(yt)
- yt=yt.T
- x=np.mat(t)# 将 X 转换为矩阵数据类型
- x=x.T# 正规方程中的 x 应该是个列向量
- w=(np.linalg.inv(x.T*x))*x.T*yt# 求解正规方程
- E_tao = np.round(1/float(w),4)# 对时间常数的 tao 的估计, 保留到 4 位小数
- E_tao
非线性函数的最小二乘拟合
说明:
1) 其间序列包含了两个等差数列 t1 和 t2 的融合, 它们的间隔互质, 没有重复, 目的是模拟采样时间的随机性. 最后用 sort()方法又对时间序列进行排序的目的是为了后续容易观察波形更直观. 如果仅仅为了使用正规方程, 是不需要重新排序的.
2) 校正的非线性方程 (3) 和原始方程 (2) 相比有一个重大的缺陷就是: 左侧求对数的括号内的数值不能为负 -- 如果是纯理论推导, 这是不可能发生的. 但假如随机噪声后的 V(t)是有可能大于阶跃幅度 Vcc 的, 此时括号内将变为一个负数, 使得计算变得没有意义. 我的解决之道是将假如的随机噪声幅度限制在仿真截止时刻 V(t)和 Vcc 之差的范围内, 及代码中的 N_amp. 由于仿真的结束时刻为 n(=5)个τ, 所以 N_amp 的值等于 np.exp(-n)*vcc.
这样做没有任何理论依据, 仅仅是受限于 (3) 和(2)式之间的非完全等价变换 -- 属于线性化校正需要付出的代价.
3) 由于待估计的参数只有一个 (1/τ) 所以正规方程的解也是只有一个元素的矩阵. 将其转换为标量后取倒得到最优估计.
三, 结果评估
加入噪声后的信号如下图所示, 与通常情况的实测波形吻合度很高.
图 2 模拟产生的带有噪声的阶跃响应
对这些加入噪声的信号进行线性校正后得到进行最小二乘估计的信号 yt 为下图所示. 其基本趋势是一条斜率为 (1/τ) 的直线, 和我预计的结果一样.
图 3 对图 2 进行线性校正后的待估计信号
但可以明显的看到, 由于 (3) 式左侧在 V(t)的大小接近 Vcc 时对加入的白噪声进行了放大. 因此当 t 递增时, 由白噪声造成的信号的不确定性大大增加了. 也就是在套用正规方程时, t 值较大时的噪声对结果的影响大于 t 值较小时的噪声对结果的影响. 这是使用非线性校正函数需要付出的重要代价.
另外, 多次运行以上代码的得到 都是一个略小于实际τ(=0.1)的数值 -- 约为 0.099 左右, 也就是说, 不是无偏估计. 这应该是由于线性校正函数 ((3) 式左侧), 对于噪声 noise 大于 0 和小于 0 的部分进行了非对称的变换造成的. 这虽然造成的偏差不大, 但也是使用非线性校正函数需要付出的代价.
四, Jupyter notebook
上述练习的一个重要目的是练习使用 Jupyter notebook, 并在其中内嵌具有很好交互性的公式等信息. 以下是部分程序运行效果的截图, 虽然我对 Markdown 语法还不熟悉, 格式和排版还不够漂亮, 但还是能够明显的提高代码的可读性.
图 4 Jupyter notebook 运行效果截图
其中需要重点记录下的是公式代码的嵌入过程:
1)我首先在纸上手写了一些公式, 用手机对其拍照, 如:
图 5 手写的公式
2)用 mathpix tools 对这些照片截图, 并扫描(mathpix tools 有 Windows 版和 iOS 版, 均可免费试用).Mathpix 可以直接得到 LaTeX 格式的公式表达式. 下图是 iOS 版本的 mathpix 界面截图.
图 6 iOS 版本的 mathpix 截图
3)在 Jupyter notebook 上部的下拉菜单中选择单元格的格式为 Markdown; 将 LaTeX 格式的公式表达式粘贴到该单元格内, 并在 LaTeX 公式表达式的前后加上 "$$" 符号, 运行该单元格即可得到图 4 所示效果的公式了.
图 7 在 Jupyter notebook 中输入 LaTeX 公式
五, 进一步的实际测试
工作做到这里离其实并没有完, 还应该做一个简单的实际电路实测一下. 我会在后续的假期中抽半天时间, 在 STM32 开发板上完成这个工作: 用 GPIO 产生一个节约信号后, 连续采集 5 个τ时间长度的信号, 并代入正规方程求解, 欢迎大家继续关注更新.
......
来源: https://www.cnblogs.com/helesheng/p/12066919.html