JP2016099104A

JP2016099104A - 温度推定システムおよび推定方法

Info

Publication number: JP2016099104A
Application number: JP2014239160A
Authority: JP
Inventors: 英朋境野; Hidetomo Sakaino
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 2014-11-26
Filing date: 2014-11-26
Publication date: 2016-05-30

Abstract

【課題】正確な室内の温度変化を推定することができるシステム等を提供すること。
【解決手段】温度推定システム１は、複数のサーバと、サーバの発熱を冷却する空調機とが設けられるデータセンタ内の温度を推定するように構成され、データセンタ内の温度推定に影響を与える実測データを取り込む入力部と、実測データに基づいて、データセンタ内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定する解法部と、モデルパラメータについて、実測データとしての温度と流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行う最適化部と、最適化されたモデルパラメータを選択して流体方程式からデータセンタ内の時間的・空間的温度変動を推定する推定部とを含む。
【選択図】図１

Description

本発明は、空間内の温度分布を推定するシステムおよび方法に関する。

サーバ室や通信機室、データセンタ等のエリアでは、空調機が設けられており、空調機がＩＣＴ(Information Communication and Technology)装置（サーバ等）の吸い込み温度を所定値以下に保つために設定温度を制御している。このような状況下において、これらのエリアでは、冷却効率を向上するため、空調機の消費電力を抑制することが重要な課題となってきている。

冷却効率を向上させるために、サーバラックの配置や冷却方法を変えて実際に実験を行うことが検討されているものの、実際の運用状態で実験するのが難しい場合もあり、数値シミュレーションによって温度分布を検証することが重要な評価手段となっている（非特許文献１、２）。

S.K. Shrivastava, J.W. VanGilder, and B.G. Sammakia, "Prediction of Cold Aisle End Airflow Boundary Conditions Using Regression Modeling," IEEE Trans. Components and Packaging Tech., vol. 30, no. 4, pp. 866-874, Dec. 2007 S.K. Shrivastava, M.I. Iyengar, B.G. Sammakia, R. Schmidt, and J.W. VanGilder, "Experimental-Numerical Comparision for a High-Density Data Center: Hot Spot Heat Fluxes in Excess of 500 W/ft2," IEEE Trans. Components and Packaging Tech., vol. 32, no. 1, pp. 166-172, March 2009

従来の数値シミュレーションでは、数理的な取り扱いの複雑さによりサーバの負荷量を一定とし、サーバ付近の温度変化だけを考慮した静的なモデルについて評価されており、時間的に変化するサーバの負荷量等の実測値については考慮されていないという問題があった。

本発明は、このような状況下においてなされたものであり、その目的は、時間的に変化する実測値を与え、この実測値を考慮することによって正確な室内の温度変化を推定するシステムおよび方法を提供することにある。

上記の課題を解決するための発明は、複数の情報処理機器と、前記情報処理機器の発熱を冷却する空調機とが設けられる対象エリア内の温度を推定する温度推定システムであって、前記対象エリア内の温度推定に影響を与える実測データを取り込む入力部と、前記実測データに基づいて、前記対象エリア内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定する解法部と、前記モデルパラメータについて、前記実測データとしての温度と前記流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行う最適化部と、前記最適化されたモデルパラメータを選択して前記流体方程式から前記対象エリア内の時間的・空間的温度変動を推定する推定部とを含む。

また、上記の課題を解決するための発明は、複数の情報処理機器と、前記情報処理機器の発熱を冷却する空調機とが設けられる対象エリア内の温度をコンピュータが推定する温度推定方法であって、前記対象エリア内の温度推定に影響を与える実測データを取り込むステップと、前記実測データに基づいて、前記対象エリア内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定するステップと、前記モデルパラメータについて、前記実測データとしての温度と前記流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行うステップと、前記最適化されたモデルパラメータを選択して前記流体方程式から前記対象エリア内の時間的・空間的温度変動を推定するステップとを含む。

本発明によれば、時間的に変化する実測値を与えることによって正確な室内の温度変化を推定することができる。

本発明の実施形態に係る温度推定システムの構成例を示すブロック図である。データセンタ内のサーバおよび空調機の配置例を示す図である。温度推定システムのシミュレーション処理例を示すフローチャートである。シミュレーションの各種設定パラメータの例を示す図である。図２の所定の上段および中段のラック上において実施されたシミュレーション結果を時系列で説明するための図である。温度推定システムで実施されるシミュレーション結果を、他のシミュレーション結果と対比して示した図である。

以下、本発明の温度推定システムの実施形態について説明する。この温度推定システム１は、データセンタ等の室内の温度分布の推定を行うものである。

[温度推定システムの構成]
図１は、本実施形態の温度推定システム１の構成例を示す図である。

図１に示すように、温度推定システム１は、ＣＰＵ（Central Processing Unit）１１と、メモリ１２と、表示部１３と、入力部１４とを含む。

ＣＰＵ１１は、各構成要素とバスで接続されて制御信号やデータの転送処理を行う。また、ＣＰＵ１１は、予測装置１全体の動作を実現するためのプログラムの実行等を行う。本実施形態のＣＰＵ１１は、下記の解法部、最適化部および推定部の機能を有するが、詳細は後記する。

解法部は、実測データに基づいて、対象エリア（データセンタなど）内の時間的・空間的温度変動を後述する熱流の流体方程式から数式的に求めるためのモデルパラメータを推定する。

最適化部は、上記モデルパラメータについて、実測データとしての温度と流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行う。

推定部は、最適化されたモデルパラメータを選択して流体方程式から対象エリア内の時間的・空間的温度変動を推定する。

メモリ１２は、ＲＯＭ(Read Only Memory)とＲＡＭ(Random Access Memory)とを含んで構成される。ＲＯＭには、温度推定装置１全体の動作制御に必要なプログラムや各種のデータ（例えば、サーバの吸込み温度、空調機の吹き出し温度などの実測データ）が記録される。ＲＡＭには、データやプログラムを一時的に記憶するための記録領域が設けられ、プログラムやデータが保持される。

表示部１３は、例えば、液晶ディスプレイ、ＥＬ(Electro-Luminescence)などのフラットパネル表示器とすることができる。

入力部１４は、後述する実測データを取り込む。入力部１４が、例えばＬＡＮ(Local Area Network)等の通信ネットワークに接続されているインターフェースである場合、入力部１４は、通信ネットワークを介して、各種センサから実測データを取り込む。

［データセンタ内の配置］
次に、この温度推定システム１によって温度分布が推定されるデータセンタ内のサーバおよび空調機の配置の一例について、図２を参照して説明する。

図２は、データセンタ内のサーバおよび空調機の配置例を示す図であって、（ａ）はデータセンタの平面図、（ｂ）はデータセンタ内部を上からみた図、（ｃ）はデータセンタ内の斜視図、（ｄ）はサーバが収納されたラック内の様子、を示す。

図２（ａ）に示すように、データセンタ内には、複数のラックが複数の列で配置される。ラックは寸法の異なるものが２つある。例えば、ラック（Ｒ３−１〜Ｒ３−７，Ｒ４−１〜Ｒ４−７）の寸法は、幅６００ｍｍ×奥行６００ｍｍであり、ラック（Ｒ５−１〜Ｒ５−６）の寸法は、幅７００ｍｍ×奥行９００ｍｍである。

また、データセンタ内には、情報処理機器としてのサーバの発熱を冷却する３台の空調機Ａ１，Ａ２，Ａ３が設けられる。

空調機Ａ１〜Ａ３は、データセンタ内を空調して室内の温度を制御する。空調機Ａ１〜Ａ３には、吹き出し温度を検知するためのセンサを備え、そのセンサ出力が温度推定装置１へ伝送されるように構成されている。空調機Ａ１〜Ａ３の吹き出し温度が検知可能であれば、センサはいかなる場所に設置し得る。なお、各空調機Ａ１〜Ａ３の寸法は、幅２１００ｍｍ×奥行９００ｍｍである。

なお、図２（ａ）の一点鎖線は境界条件となる壁を示している。図２（ｂ）および図２（ｃ）は、上記境界条件下においてデータセンタ内部のサーバおよびラックの配置を例示的に示している。図２（ｃ）において、空調機Ａ１〜Ａ３の高さは２０００ｍｍである。

各ラックは、図２（ｄ）に一例を示すように、上段、中段および下段の３段からなり、各段にそれぞれサーバ１０ａ，１０ｂ，１０ｃが収納されている。ラックの高さＨは、２１００ｍｍである。なお、図２（ｄ）では幅６００ｍｍ×奥行６００ｍｍのラックが示されているが、幅７００ｍｍ×奥行９００ｍｍのラックも３段からなり、各段にサーバが設けられるようになっている。このラックの高さも、２１００ｍｍである。なお、図２に示す様々な寸法は一例であり、変更することができることは当業者にとって自明である。

次に、このデータセンタ内の温度分布を推定するために実行される処理の全体について、図３を参照して説明する。

図３は、温度推定システム１の温度推定の処理例を示すフローチャートである。

図３において、入力部１４は、データセンタ内における実測データを取り込む（ステップＳ１００）。実測データは、例えば、サーバ吸気および排気、室内の温度、排気風量、空調機の吸い込み、排気温度、排気風量、環境温度などがある。これらの実測データの取り込みは、例えば、予め設置されている各種センサの信号を受信することによって実現される。

ステップＳ１１０では、ＣＰＵ１１は、後述する２次元または３次元の熱流体方程式の解法を計算する。この場合、数理方程式に関する条件は、事前にユーザインタフェース（UI)を通じて設定されており、ＣＰＵ１１は、この条件を使用して、上記解法を計算する。この計算の処理については、後に詳細に説明する。

なお、ステップＳ１２０において、ＣＰＵ１１は解法部として機能する。

次に、ＣＰＵ１１は、上述した実測データと数値シミュレーション値との誤差を解析して評価関数（後記式（１６））の数値を算出し（ステップＳ１２０）、モデルパラメータに関して後記遺伝的アルゴリズムを適用して最適化を実施する（ステップＳ１３０）。そして、上記最適化の判定が終了するまで（ステップＳ１４０のＹＥＳ）、ステップＳ１１０〜Ｓ１３０の最適化の判定を繰り返す。この最適化の判定では、遺伝的アルゴリズムによって複数の世代交代が繰り返され、それぞれにおけるモデルパラメータの中から、ステップＳ１２０で最適と判断されたモデルパラメータが選択される（ステップＳ１５０）。ステップＳ１５０では、ＣＰＵ１１はさらに、選択されたモデルパラメータを熱流体方程式からデータセンタ内の時間的・空間的温度変動を推定する。

なお、ステップＳ１３０〜Ｓ１４０において、ＣＰＵ１１は最適化部として機能し、ステップＳ１５０において、ＣＰＵ１１は推定部として機能する。

ＣＰＵ１１は、上記選択されたモデルパラメータを最終結果としてメモリ１２に保存した（ステップＳ１６０）後、当該データを表示部１３上に表示する（ステップＳ１７０）。

次に、上述した熱流体方程式（式（１）〜式（１５））について説明する。

先ず、熱流体方程式で使用される様々なデータに関連して、これらの表記について説明する。

以下の説明では、u，ｐ，ρ，μ，ｆは、それぞれ、速度，気圧，流体密度，粘性計数，外部力を示す。扱う対象が気流現象であるため、μをゼロとする。この場合、速度と気圧に関する支配方程式であるナビエ・ストークス（NS)における拡散項の計算を削除でき、その結果、大幅な高速解法が可能となる。

後述の方程式群は、規則的なカーテシアン格子により離散化する。気圧や速度に関する変数は、それぞれ、メッシュの中心と格子点上に配置するMACグリッドを用いる。

f_coldwind，f_hotwind，f_{thermalbuoyancy}，f_confinementは、それぞれ、空調機からの風力，サーバからの風力，熱流の浮力を示す。f_confinementは、本実施形態で採用する渦粒子法により発生する力を示す。また、T，T_amb は、それぞれ、流体の気温，環境温度を示す。

z = (0,0,1)は、３次元における鉛直方向の単位ベクトルを表記する。NS式は、修正されたMacCormackアルゴリズムにより効率がよく、精度もよい。

流体の非圧縮性を仮定して、速度の発散がゼロ（連続式）という拘束条件のもとで、NS式と連立させるべく、MICCG (0)法を投影項に適用する。ここで、ｔを時間とする。また、∇，∇²は、それぞれ空間一次微分，空間二次微分に関する演算子である。この場合、NS式は、下記式（１）で表される。

連続式は、下記式（２）で表される。

NS式における総外部力は、下記式（３）で表される。

ここで、調整係数をα，βとすると、浮力は、下記式（４）で表される。

渦粒子法による外部力は、下記式（５）で表される。

ここで、式（５）のVortex Partcle Methodは、後述する渦粒子法を表す。

気温の時間変化に関するエネルギー方程式は、下記式（６）で表される。

式（６）において、c_k は熱伝導率、T_in はサーバの吸気温度、P はサーバの負荷電力（負荷量）、をそれぞれ示す。また、C_air は空気の熱容量、Δt は時間間隔、S はサーバ吹き出し口の大きさ、V はサーバの吹き出し速度、をそれぞれ示す。なお、c_k はユーザによって定義される。T_sourceは、下記式（７）で表される。

修正されたMacCormackアルゴリズムとMICCG (0) 法は気温に関する方程式を解くときにも用いられる。

式（７）において、T_coldsource はユーザによって定義され、T_hotsource は下記式（８）で表される。

本実施形態では、サーバごとの排気温度については、サーバの熱容量，発熱量，吸気温度により決定されるものと仮定する。そして、上記定義された変数とサーバからの排気される気流の速度Vserとにより、サーバに関する気流の総質量Ｍは、下記式（９）で表される。

なお、サーバの排気温度は、上記式（８）から、下記式（１０）で与えられる。

なお、P_ser ，V_ser はともに、ユーザによって定義される。

空調機の消費電力P_crac は、下記式（１１）で与えられる。

ここで、V_acw は空調機（Ａ１〜Ａ３）からの冷気吹き出し速度、 T_{crac_in} は空調機の吸気温度，T_{crac_out} は空調機の排気温度、S_in は空調機の吹き出し部の断面積、をそれぞれ示す。V_acw，T_{crac_out} はともに、ユーザによって値が定義される。

データセンタ内の通路の風力をF _cold _Aisle、通路の温度をT _cold _Aisleとすると、これらは下記式（１２）〜（１４）で表される。

この実施形態では、実際に計測されたサーバの負荷量，吸気，排気温度，空調機の温度などに基づいて、上述した式における以下の9つのモデルパラメータを確率的なパラメータ推定法の一つである遺伝的アルゴリズム（GA)により推定する。

f_hotwind ：サーバの吹き出しによる風力
f_coldwind ：空調機の吹き出しによる風力
α , β：浮力に関するパラメータ
c_k ：熱伝導率
T_coldsource ：空調機の吹き出し冷気温度
a,b,c ：ガウス関数用(強度，平均，分散）

上記遺伝子アルゴリズム（GA）の推定処理（ステップＳ１３０）に関連して、下記ステップ１〜５の処理について説明する。

［初期値生成（ステップ１)］
初期値生成は、探索と反復の初期点を設定する。データ解の候補を遺伝子で表現した個体（N個）を用意する。各個体は未知数(f_hotwind, f_coldwind, α, β, c_k, T_coldsource, a, b, c)を表現しなおしたものとする。

GAでは現世代と次世代を想定する。まず、経験的に設定された値（現世代）に乱数を加算する。計算機で扱う上で、解をバイナリ―（０，１）列とし、これはGAでは染色体（クロモソーム）と呼ばれている。

［適合度関数処理（ステップ２)］
適合度関数処理では、適合度関数（評価関数）を用いて、現世代の各個体の適応度をそれぞれ計算する。ある確率で次の３つの動作のどれかを行い、その結果を次世代に保存する。（１）個体を２つ選択して交叉する、（２）個体を一つ選択して突然変異を行う、のうち（３）個体を一つ選択してそのままコピーする。次世代の個体数がN個になるまで上記（１）〜（３）の処理を繰り返す。そして、次世代の個体数がN個になったら次世代の内容をすべて現世代に移す。適合度関数は、下記式（１５）で定義される。

式（１５）において、分母F_mixは誤差値，Ffitnessは非負値、をそれぞれ示す。誤差が小さいほど、適合度が高くなる。即ち、現世代におけるクロモソームが選択される確率が高くなる。

［選択（ステップ３)］
この実施形態において、選択とは生物の自然淘汰を表現したものである。適合度に基づいて、個体を増やしたり消滅させたりする操作である。選択のアルゴリズムには、例えば、ルーレット選択がある。

［交叉と突然変異（ステップ４)］
交叉とは、生物が交配によって子孫を残すことをモデル化したものであり、個体の遺伝子の一部を入れ替える操作を意味する。この確率はユーザによって定義される。一点交叉、二点交叉、多点交叉が知られているが、本実施形態では、二点交叉を採用する。これは、交叉点をランダムで二つ選択し、二つの交叉点に挟まれている部分を入れ替えるものである。また、突然変異は、個体の遺伝子の一部を変化させる操作を意味する。その効果は局値に陥ることを抑制する。

［目的関数の計算（ステップ５)］
目的関数は、下記式（１６）で表される。

式（１６）において、Fmixが所定の値よりも小さくなるまで、反復的に計算が続けられる。

［渦粒子法の概要］
上述した渦粒子法は一般に、下記式（１７）〜式（２３）によって表されることが知られている（たとえば、「A Vortex Particle Method for Smoke, Water and Explosions Selle, A., Rasmussen, N. and Fedkiw, R. SIGGRAPH 2005, ACM TOG 24, p. 910-914」を参照）。

流れの中に複数の粒子を挿入し、個々の粒子には、２次元または３次元のclampedガウス分布を与えることで、局所的に流体との相互作用から擬似的な乱流表現が高速にできる。ガウス分布については、半径と強度を変化させることができる。

粒子ｐのｘでの渦力をＦｐ（ｘ）とし、格子点と粒子点とをそれぞれ、ｘとｘｐとすると、下記式が成り立つ。

ここで、式（１９）において、Ｎ_ｐ（ｘ）は計算格子点ｘの渦点の単位ベクトルを、βは調整パラメータを、それぞれ示す。Ｎ_ｐ（ｘ）は下記式（２０）で表される。

上記式（５）の渦力方程式は複数個の粒子によるものであり、下記式（２１）に置き換えることができる。なお、ガウス分布は、一定半径以上では０にclampする。

時間発展の渦度方程式は、下記式（２３）で表される。

図４は、図２に示したデータセンタ室内の温度分布を推定するために実行されるシミュレーション（図３のステップＳ１００〜Ｓ１７０の処理）時に許容できる様々なモデルパラメータの設定条件を示す図である。なお、図４は、モデルパラメータの種類を例示的に示しているに過ぎず、適宜変更することができる。

シミュレーションでは、事前に、計算格子数や計算の反復回数、収束条件などを設定することができる。

図５（ａ）および（ｂ）は、データセンタ室内の同一時点における実測値とシミュレーション結果との間の各温度の平均誤差の一例を示している。図５（ａ）に示した例では、例えば図２のラック（Ｒ３−４）上段に対する平均誤差が示され、図５（ｂ）に示した例では、例えば図２のラック（Ｒ３−４）中段に対する平均誤差が示されている。

図５（ａ）および（ｂ）の各々において、データｄ１は、ラック（Ｒ３−４）における温度の実測値が時系列に示されている。そして、データｄ２，ｄ３，ｄ４は、それぞれ、ＧＡでの反復回数を３，９，１５としたときの上記実測値とシミュレーション結果との平均誤差が示されている。

図５から、一定の時間が経つと非定常な温度変化になり、平均誤差が小さくなることがわかる。

また、上述したように、実測値とシミュレーション結果との平均誤差（温度）を求める場合において、実測値を補間することで平均誤差をさらに小さくするようにすることもできる。例えば図６は、従来の確定型最適化方法によるものと、本実施形態の確率型最適化方法（ＧＡ）によるものと、実測値を補間したものを適用したＧＡによるものとが、それぞれ例示してある。

確定型最適化法では、一般的に知られている線形二乗法に基づき、既知データ（実測値）とモデルパラメータとの関係性に基づいて室内の温度を推定することにより推定された温度を評価した。図６の例では、平均誤差が小さくなるのは、本実施形態の確率型最適化方法（ＧＡ）による場合と、実測値を補間する場合である。

一般に、元のデータ（実測値）は、データセンタ内に配置されるセンサによって検出されることになるが、このセンサの配置パターン（数、間隔、設置場所など）により、十分な実測値を得られデータがないことが多く、ＧＡによるパラメータ推定の影響を受けやすくなる。

この観点から、図６では、空間的にデータが連続的に存在することを前提として、隣接する２点間の単純な平均法やスプライン補間、ラグランジュ補間法、バイリニア法など広く用いられている手法（例えば、大宮司久明著、数値流体力学全集、第４章、補間と関数近似その１、平成26年11月13日検索、http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/Chapter4.pdf）を適用して前処理を行うのが好ましい。図６の例では、上記実測データの補間法としてバイリニア法を適用し、隣接するラックの柱上の温度から補間値を生成して与えたが、他の補間法を適用するようにしてもよい。

以上説明したように、本実施形態の温度推定システム１によると、対象エリア内の温度推定に影響を与える実測データに基づいて、対象エリア内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定し、モデルパラメータについて、実測データとしての温度と流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行う。そして、最適化されたモデルパラメータを選択して流体方程式から対象エリア内の時間的・空間的温度変動を推定する。これにより、モデルパラメータが時間変化に対応して最適化され、従って空間内の温度経時変化の推定を行うことができる。

１温度推定システム
１１ＣＰＵ
１２メモリ
１３表示部
１４入力部

Claims

複数の情報処理機器と、前記情報処理機器の発熱を冷却する空調機とが設けられる対象エリア内の温度を推定する温度推定システムであって、
前記対象エリア内の温度推定に影響を与える実測データを取り込む入力部と、
前記実測データに基づいて、前記対象エリア内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定する解法部と、
前記モデルパラメータについて、前記実測データとしての温度と前記流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行う最適化部と、
前記最適化されたモデルパラメータを選択して前記流体方程式から前記対象エリア内の時間的・空間的温度変動を推定する推定部と
を含むことを特徴とする温度推定システム。
前記実測データは、前記情報処理機器の負荷量、吸気温度、排気温度、および前記対象エリア内の環境温度を含むことを特徴とする請求項１に記載の温度推定システム。
前記流体方程式は、ナビエ・ストークス方程式を含むことを特徴とする請求項１に記載の温度推定システム。
複数の情報処理機器と、前記情報処理機器の発熱を冷却する空調機とが設けられる対象エリア内の温度をコンピュータが推定する温度推定方法であって、
前記対象エリア内の温度推定に影響を与える実測データを取り込むステップと、
前記実測データに基づいて、前記対象エリア内の時間的・空間的温度変動を熱流の流体方程式から数式的に求めるためのモデルパラメータを推定するステップと、
前記モデルパラメータについて、前記実測データとしての温度と前記流体方程式から数式的に求められる温度との誤差が所定値以下となるよう、遺伝子アルゴリズムによる最適化を行うステップと、
前記最適化されたモデルパラメータを選択して前記流体方程式から前記対象エリア内の時間的・空間的温度変動を推定するステップと
を含むことを特徴とする温度推定方法。
前記実測データは、前記情報処理機器の負荷量、吸気温度、排気温度、および前記対象エリア内の環境温度を含むことを特徴とする請求項４に記載の温度推定方法。
前記流体方程式は、ナビエ・ストークス方程式を含むことを特徴とする請求項４に記載の温度推定方法。