常用算法之二牛顿迭代法求解PT100温度（高阶方程求解）

雨点打透心脏的1/2处 2022-04-16 05:15 695阅读 0赞

# 问题 #

在最近的工作中用到了PT100铂电阻，它适用于医疗、电机、工业、温度计算、卫星、气象、阻值计算等高精温度设备，应用范围非常之广泛。类似的还有PT1000。  
  PT后的100即表示它在0℃时阻值为100欧姆，在100℃时它的阻值约为138.5欧姆。它的工作原理：当PT100在0摄氏度的时候他的阻值为100欧姆，它的的阻值会随着温度上升它的阻值是成匀速增长的。  
  在国家标准《GBT 30121-2013 工业铂热电阻及铂感温元件.pdf》中（下载见附件），PT100的分度表如下图所示：  
![PT100][]  
温度与阻值的计算公式，在标准文档中也是有的，如下：  
![PT100Func][]  
以上公式涉及到高阶方程的求解问题！

# 方法 #

一般可使用牛顿迭代法求近似解！关于该方法，网友[马同学高等数学 ][Link 1]写了一篇非常棒的文章[如何通俗易懂地讲解牛顿迭代法？][Link 2]来进行了非常详细的介绍。由于没有作者的授权，这里不方便转载，感兴趣的可以点击以上链接去看看！*需要很强的数学知识哦！！！突然感觉，当初的高等数学白学了，唉。。。*

# 源码 #

从温度获取阻值非常简单，重点在于从阻值获取温度！用二阶导数切线法，根据电阻接近零度电阻程度，用1到3次计算可得到较为精确的温度值：

1.  先用线性公式得到一个粗略的温度：  
    t1 = (Rt / R0 - 1) / A  
    计算t1的阻值：  
    0度以上用公式：Rt1 = R0 \* (1 + A \* t1 + B \* t1 \* t1);  
    0度以下用公式：Rt1 = R0 \* \[1 + A \* t1 + B \* t1 \* t1 + C \* (t1 - 100) \* t1 \* t1 \* t1\];  
    如果 |Rt1 - Rt| < 0.001，t1 就是所要的温度，反之继续进行下面的计算：
2.  计算一阶导数和二阶导数：  
    如果 Rt ≥ R0  
    t1’ = 1 / \[R0 \* (A + 2 \* B \* t1)\]  
    t1’’ =-2 \* B \* R0 \* t1’ \* t1’ \* t1’  
    如果 Rt < R0  
    t1’ = 1 / \[R0 \* (A + 2 \* B \* t1 - 300 \* C \* t1 \* t1 + 4 \* C \* t1 \* t1 \* t1)\]  
    t1’’=- R0 \* (2 \* B - 600 \* C \* t1 + 12 \* C \* t1 \* t1) \* t1’ \* t1’ \* t1’
3.  基于 Rt，t1，Rt1 计算近似温度 t2：  
    t2 = t1 + t1’ \* (Rt - Rt1) + 0.5 \* t1’’ \* (Rt - Rt1) \* (Rt - Rt1)，再计算出 t2 对应的电阻 Rt2。
4.  如果 |Rt2 - Rt| < 0.001，t2 就是所要的温度，反之再从第二步开始计算 t2 的一二阶导数，最终得出精确的温度值。

#include <stdio.h>
    #include <math.h>
    #include <stdbool.h>
    
    #define A 3.9083e-3
    #define B -5.775e-7
    #define C -4.183e-12
    
    
    /** * @brief 从温度获取阻值 * @param[in] temp 温度值 * @param[out] resist 阻值 * @retval 是否成功 */
    bool GetResistanceFromTempereture(double temp, double* resist)
    { 
    	double fT, fR;
    	bool ret;
    
    	if(resist == NULL)
    	{ 
    		return false;
    	}
    
    	if (temp >= -200 && temp < 0)
    	{ 
    		*resist = 100 * (1 + A * temp + B * temp * temp + C * (temp - 100) * temp * temp * temp);
    		ret = true;
    	}
    	else if (temp >= 0 && temp <= 850)
    	{ 
    		*resist = 100 * (1 + A * temp + B * temp * temp);
    		ret = true;
    	}
    	else
    	{ 
    		*resist = 0;
    		ret = false;
    	}
    	return ret;
    }
    
    /** * @brief 从阻值获取温度 * @param[in] resist 阻值 * @param[out] temp 温度值 * @retval 是否成功 * @note 用二阶导数切线法，根据电阻接近零度电阻程度，用1到3次计算可得到较为精确的温度值： * 1. 先用线性公式得到一个粗略的温度： * t1 = (Rt / R0 - 1) / A * * 计算t1的阻值： * 0度以上用公式：Rt1 = R0 * (1 + A * t1 + B * t1 * t1); * 0度以下用公式：Rt1 = R0 * [1 + A * t1 + B * t1 * t1 + C * (t1 - 100) * t1 * t1 * t1]; * * 如果 |Rt1 - Rt| < 0.001，t1 就是所要的温度，反之继续进行下面的计算： * * 2. 计算一阶导数和二阶导数： * 如果 Rt ≥ R0 * t1' = 1 / [R0 * (A + 2 * B * t1)] * t1'' =-2 * B * R0 * t1' * t1' * t1' * 如果 Rt < R0 * t1' = 1 / [R0 * (A + 2 * B * t1 - 300 * C * t1 * t1 + 4 * C * t1 * t1 * t1)] * t1''=- R0 * (2 * B - 600 * C * t1 + 12 * C * t1 * t1) * t1' * t1' * t1' * * 3. 基于 Rt，t1，Rt1 计算近似温度 t2： * t2 = t1 + t1' * (Rt - Rt1) + 0.5 * t1'' * (Rt - Rt1) * (Rt - Rt1)，再计算出 t2 对应的电阻 Rt2。 * * 4. 如果 |Rt2 - Rt| < 0.001，t2 就是所要的温度，反之再从第二步开始计算 t2 的一二阶导数，最终得出精确的温度值。 */
    bool GetTemperetureFromResistance(double resist, double* temp)
    { 
    	bool ret;
    	double fT, fT0;
    	short i;
    
    	if(temp == NULL)
    	{ 
    		return false;
    	}
    
    	/* 1. 先用线性公式得到一个粗略的温度 */
    	fT0 = (resist / 100 - 1) / A;
    	/* 计算t1的阻值 */
    	if (resist >= 18.52 && resist < 100) /* -200℃ ~ 0℃ */
    	{ 
    		for (i = 0; i < 50; i++)
    		{ 
    			fT = fT0 + (resist - 100 * (1 + A * fT0 + B * fT0 * fT0 - 100 * C * fT0 * fT0 * fT0 + C * fT0 * fT0 * fT0 * fT0)) /
    						   (100 * (A + 2 * B * fT0 - 300 * C * fT0 * fT0 + 4 * C * fT0 * fT0 * fT0));
    			if (fabs(fT - fT0) < 0.001) /* 如果 |Rt1 - Rt| < 0.001，t1 就是所要的温度 */
    			{ 
    				break;
    			}
    			else
    			{ 
    				fT0 = fT;
    			}
    		}
    		*temp = fT;
    		ret = true;
    	}
    	else if (resist >= 100 && resist <= 390.481) /* 0℃ ~ 850℃ */
    	{ 
    		for (i = 0; i < 50; i++)
    		{ 
    			fT = fT0 + (resist - 100 * (1 + A * fT0 + B * fT0 * fT0)) / (100 * (A + 2 * B * fT0));
    			if (fabs(fT - fT0) < 0.001) /* 如果 |Rt1 - Rt| < 0.001，t1 就是所要的温度 */
    			{ 
    				break;
    			}
    			else
    			{ 
    				fT0 = fT;
    			}
    		}
    		*temp = fT;
    		ret = true;
    	}
    	else
    	{ 
    		ret = false;
    	}
    	return ret;
    }
    
    /* 测试 */
    int main(int argc, char const *argv[])
    { 
    	double temp = 0;
    	double resist = 68.33;
    
    	if(GetTemperetureFromResistance(resist, &temp))
    	{ 
    		printf("Temp = %f\r\n", temp);
    	}
    	else
    	{ 
    		printf("Temp = Error!\r\n");
    	}
    	return 0;
    }

# 附件 #

1.  [GBT 30121-2013 工业铂热电阻及铂感温元件.pdf][GBT 30121-2013 _.pdf]

[PT100]: /images/20220416/8201d72de6dd4dcaa5951b31840d7867.png
[PT100Func]: /images/20220416/ac13917aa06a436fb54e284ec6a708bb.png
[Link 1]: https://www.matongxue.com
[Link 2]: https://www.matongxue.com/madocs/205.html
[GBT 30121-2013 _.pdf]: https://download.csdn.net/download/zcshoucsdn/10790573