基于MATLAB的冷却液温度传感器特性曲线拟合方法研究

Research on Fitting Method of Characteristic Curve of Coolant Temperature Sensor Based on MATLAB

摘要: 本文通过对比各曲线拟合方法的优缺点，选择最优的拟合方法：最小二乘法，并利用MATLAB软件中集成的最小二乘法，对比各阶次最小二乘法特性曲线拟合结果，得到较理想的冷却液温度传感器特性曲线方程。对比结果表明五阶次拟合得到的特性曲线效果较理想。

Abstract: The best fitting method: the least square method is selected by comparing the advantages and dis-advantages of curve fitting in this paper and the least square method integrated in MATLAB soft-ware is used to compare the results of characteristic curve fitting of each order of least square method, so that the characteristic curve equation of coolant temperature sensor is closer to the real one. The results show that the effect of the characteristic curve obtained by the fifth order fitting is ideal.

1. 引言

在整车故障模拟与消除时，有时会通过信号发生器发送电压信号来模拟传感器反馈给整车控制单元的信号，其中最常见的就是通过模拟冷却液温度信号，来实现冷却液温度相关故障的产生与清除。如果得到冷却液传感器特性曲线方程，就可以通过发送电压信号精确控制控制单元识别出的冷却液温度。得到特性曲线方程的关键是选择合适的曲线拟合方法。王荣忠等利用单片机提供算力，采用最小二乘法进行曲线拟合 [1] 。因最小二乘法涉及复杂数学推导过程以及复杂数学公式，上述方法需大量C语言编程工作，实现较为复杂，工作量大。本文首先介绍曲线拟合基本知识，确定拟合用到的方法，结合MATLAB自带数学工具，通过不同阶次拟合结果对比，得出最优拟合方式，研究方法简单、稳定、可靠，为基于MATLAB的冷却液特性曲线拟合方法提供思路和方法。

2. 数据拟合理论基础

2.1. 曲线拟合数值理论

曲线拟合(curve fitting)是一种数据处理方法，其用不间断的曲线来反映坐标系中的一些给定的离散点的总体趋势，这些坐标点间存在一定的函数关系。在科学计算中，由实验得到x与y的数据集 $\left(x_i,y_i\right)$ 。人们希望通过使用一种与离散坐标点的规律相一致的解析式来表示x与y之间的对应函数关系，即在一定的条件下最优地逼近已给定的离散数据坐标点 [2] 。

由于存在的噪声，在拟合曲线时，拟合得到的曲线并不一定要经过每个实验给定的数据点，而只要求拟合出来的函数与数据点的总体误差最小就可以了。通常将近似的值 $s\left(x_i\right)$ 与实验测定值 $y_i$ 作差，差值 ${\delta }_i$ 就是所谓的残余误差，即 ${\delta }_i=s\left(x_i\right)-y_i$ 。一定程度上，残余误差结果的大小被视为判断拟合运算成功与否的主要标志之一。目前公认的三个可以衡量拟合效果的准则为：1) 最大残余误差对应的绝对值应该最小： $\max \left|{\delta }_i\right|=\min$ ；2) 残余误差对应的各绝对值的和应该是最小的： $\sum \left|{\delta }_i\right|=\min$ ；3) 残余误差对应的平方和应该是最小的： $\sum {\delta }_i^2=\min$ 。对于前述的三个准则，准则(1)与(2)是不易在实际中计算的，因为两者都存在绝对值，因此我们更多的采纳第三条准则。通过第三条准则来计算某些参数，从而获得拟合所需曲线的算法就是最小二乘法。即最小二乘法是一种数值优化技术，通过使误差平方和最小化，从而得到最佳匹配函数 [3] 。其计算相对简单实用，结果较为准确，是目前实验中常用的函数逼近方法。

2.2. 最小二乘法数值算法

假设通过实验给定的离散数据点为 $\left\{\left(x_i,y_i\right),i=1,2\cdots ,m\right\}$ ，数据点对应的权系数 $a_i$ ，从而得到自变量x和因变量y对应的函数关系 $y=s\left(x;a_0,a_1,\cdots ,a_n\right)\left(n<m\right)$ ，最小二乘法并不一定要 $y=s\left(x\right)$ 经过每一个给定数据点，但要保证对应的残差最小。设逼近函数为：

$s\left(x\right)=a_0{\phi }_0\left(x\right)+a_1{\phi }_1\left(x\right)+\cdots +a_n{\phi }_n\left(x\right)$ (1)

假定已有一些数据点 $\left\{\left(x_i,y_i\right),i=0,1,2,\cdots ,m\right\}$ 以及对应的一组权：

$\left\{{\rho }_i,i=0,1,2,\cdots ,m\right\}$

求 $y=s\left(x;a_0,a_1,\cdots ,a_n\right)\left(n<m\right)$ 使

${\Vert \delta \Vert }_2^2=p\left(a_0,a_1,\cdots ,a_n\right)={\displaystyle \underset{i=0}{\overset{m}{\sum }}{\rho }_i}{\left[s\left(x_i\right)-y_i\right]}^2$ (2)

的值最小，则 $\left(a_0^{\ast },a_1^{\ast },\cdots ,a_n^{\ast }\right)$ 是多元函数 $p\left(a_0^{},a_1^{},\cdots ,a_n^{}\right)$ 得到极小值时对应的点，这就是所谓的最小二乘法 [4] 。由求取极值所需的一些必要条件能够得到：

$\frac{\partial p}{\partial a_k}=2{\displaystyle \underset{i=0}{\overset{m}{\sum }}{\rho }_i}\left[a_0{\phi }_0\left(x_i\right)+a_1{\phi }_1\left(x_i\right)+\cdots +a_n{\phi }_n\left(x_i\right)-y_i\right]{\phi }_k\left(x_i\right),k=0,1,\cdots ,n$ (3)

依据内积的定义可以引入对应的带权内积记号：

$\left\{\begin{array}{l}\left({\phi }_j,{\phi }_k\right)={\displaystyle \underset{i=0}{\overset{m}{\sum }}{\rho }_i{\phi }_j\left(x_i\right){\phi }_k\left(x_j\right)}\\ \left(y,{\phi }_k\right)={\displaystyle \underset{i=0}{\overset{m}{\sum }}{\rho }_i{y}_i{\phi }_k\left(x_j\right)}\end{array}\right\}$ (4)

则式(3)可改写为：

$\left({\phi }_{0,}{\phi }_k\right)a_0+\left({\phi }_{1,}{\phi }_k\right)a_1+\cdots +\left({\phi }_n{\phi }_k\right)a_n=\left(y,{\phi }_k\right),k=0,1,\cdots ,n$ (5)

这是关于参数 $a_0,a_1,\cdots ,a_n$ 的线性方程组，用矩阵的形式表示为：

$\left[\begin{array}{l}\left({\phi }_0,{\phi }_0\right)\left({\phi }_0,{\phi }_0\right)\cdots \left({\phi }_0,{\phi }_0\right)\\ \left({\phi }_0,{\phi }_0\right)\left({\phi }_0,{\phi }_0\right)\cdots \left({\phi }_0,{\phi }_0\right)\\ ⋮⋮\cdots ⋮\\ \left({\phi }_0,{\phi }_0\right)\left({\phi }_0,{\phi }_0\right)\cdots \left({\phi }_0,{\phi }_0\right)\end{array}\right]$ (6)

由于 ${\phi }_0\left(x\right),{\phi }_1\left(x\right),\cdots ,{\phi }_n\left(x\right)$ 线性无关，故方程组存在唯一的解 $a_k=a_k^{\ast },k=0,1,\cdots ,n$ 于是有：

$\stackrel{*}{p}\left(x\right)={\displaystyle \underset{i=0}{\overset{m}{\sum }}{a}_k}\stackrel{\ast }{{\phi }_k}\left(x\right)$ (7)

3. MATLAB中的最小二乘法

MATLAB是Math Works公司推出的一种集数值计算、符号运算及出色的图形处理、程序语言设计等功能于一体的科学计算软件 [5] 。

前章所述算法即为最小二乘法的数值计算的理论基础，本文通过在MATLAB中使用最小二乘法对冷却液传感器特性曲线进行拟合，确定最优的拟合阶次，从而得到特性曲线的方程。用MATLAB对给定的数据点作处理，只需要利用软件内部集成的程序，编制很简单的程序文本，运行并编译程序，就能够获得所需要的结果。这样就解决了最小二乘法的很多缺点，比如计算量大、计算枯燥等问题，同时使得繁琐的数学运算转换为简单的操作，并且能够精确地绘制出所需要的拟合曲线 [6] [7] 。

下文简单介绍通过MATLAB实现最小二乘法处理实验数的方法。最小二乘的算法是：

1) 导入实验数据点 $\left(x_i,y_i\right)$ 与数据点对应的权系数 $\left\{{\rho }_i,i=0,1,2,\cdots ,m\right\}$ ，选取逼近 ${\left\{{\phi }_k\left(x\right)\right\}}_{k=0}^{\infty }$ ；

2) 求取方程组6；

3) 由公式7可以获得函数所需最优的逼近。

将上面所叙述算法编译成为MATLAB可以运行的程序，计算得到数据点最小二乘对应的拟合系数， $\text{FunctionS}=\text{squar\_least}\left(x,y,n,w\right)$ 。

数据点的最小二乘法拟合，其中x，y为实验数据点的 $\left(x,y\right)$ 坐标；n为数据点拟合所对应的次数，缺省时 $n=1$ ；w为对应的权系数，缺省时对应的值也是1；s是给定的实验数据点作拟合运算时的系数。

4. 样本数据点曲线拟合与分析

4.1. 各阶次曲线拟合

为了进行曲线拟合，我们需要测量一组冷却液温度与电压的基础样本数据。温度(℃) {−30, −20, −10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150}，电压(V) {4.59, 4.389, 4.064, 3.657, 3.168, 2.659, 2.203, 1.76, 1.371, 1.076, 0.831, 0.642, 0.498, 0.393, 0.305, 0.241, 0.192, 0.155, 0.083}。

利用上述样本数据点，通过MATLAB软件拟合，分别得到各阶次曲线拟合的结果，具体如下：

图1为冷却液温度与电压一阶拟合曲线，拟合得到的方程如下：

y = p1x + p2 (p1 = −0.0268, p2 = 3.313)，一阶拟合的SSE (和方差)为3.976。

Figure 1. First-order temperature-voltage fitting curve

图1. 一阶温度–电压拟合曲线

图2为冷却液温度与电压二阶拟合曲线，拟合得到的方程如下：

y = p1x² + p2x + p3 (p1 = 1.634e−04, p2 = −0.0465, p3 = 3.411)，二阶拟合的SSE (和方差)为0.3538。

Figure 2. Second-order temperature-voltage fitting curve

图2. 二阶温度–电压拟合曲线

图3为冷却液温度与电压三阶拟合曲线，拟合得到的方程如下：

y = p1x³ + p2x² + p3x + p4 (p1 = 5.539e−07, p2 = 6.37e−05, p3 = −0.0435, p4 = 3.47)，三阶拟合的SSE (和方差)为0.2596。

Figure 3. Third-order temperature-voltage fitting curve

图3. 三阶温度–电压拟合曲线

图4为冷却液温度与电压四阶拟合曲线，拟合得到的方程如下：

y = p1x⁴ + p2x³ + p3x² + p4x + p5 (p1 = −1.908e−08, p2 = 5.134e−06, p3 = 2.026e−04, p4 = −0.04452, p5 = 3.618)，四阶拟合的SSE (和方差)为0.01476。

Figure 4. Fourth-order temperature-voltage fitting curve

图4. 四阶温度–电压拟合曲线

图5为冷却液温度与电压五阶拟合曲线，拟合得到的方程如下：

y = p1x⁵ + p2x⁴ + p3x³ + p4x² + p5x + p6 (p1 = 9.434e−11, p2 = −4.739e−08, p3 = 7.602e−06, p4 = 2.394e−04, p5 = 3.64)，五阶拟合的SSE (和方差)为0.002601。

Figure 5. Fifth-order temperature-voltage fitting curve

图5. 五阶温度–电压拟合曲线

4.2. 曲线拟合结果分析

对比4.1章节中各阶次曲线拟合的结果可得出：一阶、二阶、三阶拟合曲线明显未能通过大部分样本数据点，三者SSE均比较高，波形失真严重，四阶、五阶拟合曲线基本均能通过样本数据点，能够反映样本数据点的真实情况，但五阶SSE更小，即对于冷却液传感器特性曲线的拟合，五阶拟合的效果相对较好(六阶以上的拟合效果基本与五阶相同，本文不再列出)。即我们可以通过多项式y = p1x⁵ + p2x⁴ + p3x³ + p4x² + p5x + p6 (p1 = 9.434e−11, p2 = −4.739e−08, p3 = 7.602e−06, p4 = 2.394e−04, p5 = 3.64)，输入想要的冷却液温度，相对精确的计算出所对应的电压信号，更好的控制冷却液温度，为冷却液温度相关故障的模拟与消除提供一种有效便捷的电压计算方法。

5. 结论

本文使用最小二乘法，利用MATLAB软件十分便捷的最小二乘法交互式程序，得出最优的冷却液特性曲线方程，并确定五阶次拟合得到的特性曲线效果较理想。由此可见运用普通最小二乘法原理能从实验中得出的一大堆看上去杂乱无章的数据中找出一定规律，拟合成一条直线或曲线来反映所给数据点总趋势，提供了一种非常方便实效的数据处理方法。

当然，需要注意的是本文中的样本数据来源于某一型号冷却液传感器，其他种类的冷却液温度传感器得到的基础样本数据可能与本文中的不同，即传感器的特性曲线可能不同，但拟合方式与拟合阶次的选择的方法是相同的。