Arduinoで体温計をつくるその2

その1で温度の測定まではできた。
ここから温度の予測を行ってみようと思う。

Arduinoで体温測定 - YouTube
youtu.be

結論から言うとだいぶ微妙な感じになってしまった。

理論式から予測温度を導出する。

理論式をいろいろ調べてみるとニュートンの冷却の法則が近いようだ。

[http://topicmaps.u-gakugei.ac.jp/phys/matsuura/lecture/heat/contents/heatcdt.htm 熱伝導による冷却

この式と微分の関係式から体温の予測をしようと思う

温度:y[℃] と時間t[秒]の関係式が　　y=(T0-A)*exp(B*t)+A
（t=0の時の温度をT0[℃]、Aは∞秒後の温度(予測温度)、Bは表面積や比熱によって定まる定数であると仮定する）
だとすると、

y=(T0-A)*exp(B*t)+A
${\frac{dy}{dt}}$ =(T0-A)*B*exp(B*t)　　　　…①
${\frac{d^2y}{dt^2}}$ =(T0-A)*B^2*exp(B*t)　…②

②式に①式を代入して
${\frac{d^2y}{dt^2}}$ = ${\frac{dy}{dt}}$ *B…③
と表せる。

時間を横軸、温度を縦軸にとると、
${\frac{dy}{dt}}$ は時間と温度の傾きから求められるので、
${\frac{dy}{dt}}$ = ${\frac{温度の増加量}{時間の増加量}}$ …④
同様に2階微分は1階微分の傾きから得られるので
${\frac{d^2y}{dt^2}}$ = ${\frac{温度の傾きの増加量}{時間の増加量}}$ …⑤

③④⑤から定数Bを求めることができ、
Bの値と、1組の温度と時間のデータをもとの式に代入することでA(予測温度)を取得できる。

ソースコード

本来は実測温度がほぼ一定になる前に予測値を出せるはずでしたが、誤差が大きいようで精度が悪いです。温度の上昇が少なくなればそこそこな値が出ますが、そのころには測定温度がほぼ体温の値になっているので、意味がない。いろいろ考えて式の導出をやったので今の私の限界・・・

#include <math.h>

int analog = 0;          //アナログpin　生値
int V_0 = 3276;          //回路電源[mV]
double V_R0 =1;          //部分電圧[mV]
int circuit_R0 =10000;  //回路内抵抗[Ω]

int thermo_B = 3380;      //サーミスタB定数
int thermo_R0 = 10000;   //サーミスタ基準抵抗値[Ω]
int thermo_T0 = 298;     //サーミスタ基準温度[K]
double thermo_R = 1;     //サーミスタ測定抵抗値[Ω]
double temp = 25.00;     //測定温度[℃]

double temp_0=1.0;     //温度初期値
int time_0=0;          //時間初期値

int i=0;               //測定回数
double temp_data[3];  //測定温度
double time_data[3];  //経過時間

double dy_1=0.0;      //二階微分導出用一階微分の値1
double dy_2=0.0;      //二階微分導出用一階微分の値2
double dy=0.0;        //dy/dt
double d2y=0.0;       //(d/dt)*(dy/dt)二階微分の値
double C2=0.0;

double temp_forecast; //予測温度
int data=0;
double new_forecast;

void setup () {
  Serial.begin(9600);          /* 9600 bpsで接続 */
  Serial.println("READ START!");   //温度の表示
  analog=analogRead(0); 
  time_0 = micros(); 
  V_R0 = analog*3.3/ 1.024;                    //回路内抵抗の消費電圧を算出
  thermo_R=V_0/V_R0* circuit_R0 - circuit_R0;  //サーミスタの抵抗値を算出
  temp_0=(1000/(1/(0.001*thermo_T0)+log(thermo_R/thermo_R0)*1000/thermo_B)-273);  //温度の計算
}

void loop () { 
  if(i<30){
    analog=analogRead(0);                        //3回測定した電圧の平均値から温度を求める
   // analog+=analogRead(0);
    time_data[i%3] = (micros()-time_0);     
   // analog+=analogRead(0);
   // analog/=3;  
    V_R0 = analog*3.3/ 1.024;                    
    thermo_R=V_0/V_R0* circuit_R0 - circuit_R0; 
    temp_data[i%3]=(1000/(1/(0.001*thermo_T0)+log(thermo_R/thermo_R0)*1000/thermo_B)-273); 
    Serial.print("now_temperature:\t");
    
    if(i%3==2){                                  //3回温度を求めたらそのデータから予測温度を求める
      dy_1=(temp_data[1]-temp_data[0]+1.0e-80)*1000000/(time_data[1]-time_data[0]);
      dy_2=(temp_data[2]-temp_data[1]+1.0e-80)*1000000/(time_data[2]-time_data[1]);
      dy=(temp_data[2]-temp_data[0]+1.0e-80)*1000000/(time_data[2]-time_data[0]);
      d2y=2*(dy_2-dy_1)*1000000/(time_data[2]-time_data[0]);
      C2=d2y/dy;
      temp_forecast=(temp_data[1]-temp_0*exp(C2*time_data[1]))/(1-exp(C2*time_data[1]));
      if(temp_forecast>30){        //直近の算出できた予測温度を取りためる
        Serial.print(temp_data[i%3],1);   //測定温度の表示 
        Serial.print("\tforecast:\t");
        Serial.println(temp_forecast,1); 
        new_forecast=temp_forecast;
        data++;  
      }else{
        Serial.println(temp_data[i%3],1);   
     }
      delay(400);
    }else{
        Serial.println(temp_data[i%3],1);
     }
    i++;
  delay(500);
  }else{
   Serial.println("RESULT");   //最終的な予測温度の表示  
   Serial.print("forecast:\t");
   Serial.println(new_forecast,1);  
   Serial.print("now temperature:\t");
   analog=analogRead(0); 
   V_R0 = analog*3.3/ 1.024;                    
   thermo_R=V_0/V_R0* circuit_R0 - circuit_R0; 
   temp_data[2]=(1000/(1/(0.001*thermo_T0)+log(thermo_R/thermo_R0)*1000/thermo_B)-273); 
   Serial.println(temp_data[2],1);
   delay(10000);
 }
}

問題点

色々なところで誤差が出ていて、精度がよろしくない。考えられる原因をどんどん挙げてく

温度を求めるまでの誤差

明らかに測定値にばらつきがあり、誤差はあるがそこまで大きな誤差ではないようだ。
サーミスタは温度によって抵抗値が変わる。これを電圧の値から逆算して温度を求めている。
温度の値を求める段階までで起こる誤差として、電圧自体の測定誤差と、電圧から温度を求めるときのまるめ誤差などの計算による伝搬こり徐々に大きな誤差になっていることが考えられる。
他にも電流を流すことでサーミスタ自身の発熱も起こってしまうので厳密にはここでも誤差がでているだろう

予測値の導出までの誤差

求めた温度が絶対に正しいとして導出しているので、ここでも誤差の伝搬によってさらに値がずれてくる。
また理論値では、温度は単調に増加(又は減少)するが、微小時間が経過しても温度の測定値は同じ値になってしまうことなどから計算途中に0で割ることなどが起こってしまう
それを回避するために理論式に数値を加えてしまっているので、そもそもの温度が正しく計測できていても誤差が出てしまう。
さらに温度の測定と時間の計測が微妙にずれている、おそらく一番の原因である、傾きの大きさを増加量から導出するとずれが生じる、などが積もりに積もってうまく値が出ないのだと考えられる

理論値そのものの誤差

定数Bは、測温部と接触している表面積の時間変化などには対応していない。短い時間ごとにBの値を求め直しているが、それでもこの部分にも誤差が出ている。

最小二乗法とかを使ってみようとしたがなかなかうまいこといかなかったので、ひとまずこれでおわりにして、またしばらくたったら挑戦してみようと思う