"温室に配置される暖房機のONOFF制御を行う暖房機の制御方法において、制御タイミングとなる各ステップでのON又はOFFの組み合わせである複数の運転パターンについて、前記温室の物理モデルに基づいて前記温室の温度変化を予測し、予測された前記温度変化に基づいて、前記温室の温度が目標温度に近づくように、前記複数の運転パターンの中から最適運転パターンを選択し、前記最適運転パターンに基づいて前記暖房機の前記ONOFF制御を行う、ことを特徴とする暖房機の制御方法。予め定められた直近の予測区間の前記複数の運転パターンについて前記温室の温度変化を予測することを特徴とする請求項1記載の暖房機の制御方法。前記ステップごとに又は複数の前記ステップごとに前記複数の運転パターンについて前記温度変化を予測し、前記最適運転パターンのうち前記温度変化の予測間隔分に基づいて前記ONOFF制御を行う、ことを特徴とする請求項1又は請求項2記載の暖房機の制御方法。前記暖房機の制約条件を満たす複数の運転パターンについて前記温室の温度変化を予測することを特徴とする請求項1から請求項3のいずれか1項記載の暖房機の制御方法。使用頻度に応じて限定された複数の運転パターンについて前記温室の温度変化を予測することを特徴とする請求項1から請求項4のいずれか1項記載の暖房機の制御方法。前記暖房機は、木質燃料を燃焼させることで前記温室を暖房することを特徴とする請求項1から請求項5のいずれか1項記載の暖房機の制御方法。予測された前記温度変化の誤差と予め定められた燃費の重みとを含む評価関数に基づいて、前記温室の温度が目標温度に近づくように且つ燃費の悪化を抑えるように、前記複数の運転パターンの中から前記最適運転パターンを選択することを特徴とする請求項1から請求項6のいずれか1項記載の暖房機の制御方法。温室に配置され、制御部によりONOFF制御される暖房機において、前記制御部は、制御タイミングとなる各ステップでのON又はOFFの組み合わせである複数の運転パターンについて、前記温室の物理モデルに基づいて前記温室の温度変化を予測し、予測された前記温度変化に基づいて、前記温室の温度が目標温度に近づくように、前記複数の運転パターンの中から最適運転パターンを選択し、前記最適運転パターンに基づいて前記暖房機の前記ONOFF制御を行う、ことを特徴とする暖房機。上記のA~Cの特徴は、補償器設計を難しくしている原因であり、特に上記Cの特徴である運転制約は、高精度な温度制御には大きなネックとなっていた。また、上記Bの特徴である時定数が大きい点は、暖房機の時定数が温室の温度特性のそれにかなり近い大きさとなっており、非常に制御しにくい制御対象となっている。これらの特徴を有する木質ペレット焚き暖房機に対し、従来は、与えられた温室の目標温度の上側及び下側に、それぞれ暖房機消火のためのOFF点温度及び暖房機点火のためのON点温度を設定し、温室温度が設定温度に到達した時点で暖房機への運転信号を決定する方法が採られている。さらに、稼動中に運転指令を送った直後の温室温度のオーバシュート量及びアンダシュート量に応じてON点温度及びOFF点温度を微小な範囲で随時更新する手法を採用することで制御精度が改善されている。しかし、このような暖房機には以下の問題点が存在する。1与えられた目標温度に対し、ON点温度及びOFF点温度の設定基準が不明確である。2暖房機の立上り及び立下り時間時定数の情報が活かされていない。3暖房機の運転制約の影響が考慮されていない。4外気温の情報が考慮されていない。特に、温室から奪われる熱量は外気温度と温室温度との差に依存するため、ONOFF点温度の更新を併用した場合であっても対症療法的な制御であり、また本質的に室温に対する暖房機と外気温の動特性、及び制約条件を考慮していない点で制御性能は満足できるものではない。実際に、図25に示すような従来の木質ペレット焚き暖房機による温度制御実績は、目標温度である18°C周りで5°Cppピークピーク値程度であり、受粉時期における2°Cpp以内といったレベルの温度精度要求には全く対応できていない状況にある。なお、図25に示す例では、暖房機は、例えば19時~翌日8時など、外気温の下がる夜間に稼働するように制御されている。なお、上述の事項のうち特に1~4については、木質ペレット焚き暖房機に限らず、ONOFF制御される暖房機であれば他の燃料を用いたものについても同様のことがいえる。本発明の目的は、暖房機のONOFF制御を高精度に行うことができる暖房機の制御方法及び暖房機を提供することである。本発明の暖房機の制御方法は、温室に配置される暖房機のONOFF制御を行う暖房機の制御方法において、制御タイミングとなる各ステップでのON又はOFFの組み合わせである複数の運転パターンについて、上記温室の物理モデルに基づいて上記温室の温度変化を予測し、予測された上記温度変化に基づいて、上記温室の温度が目標温度に近づくように、上記複数の運転パターンの中から最適運転パターンを選択し、上記最適運転パターンに基づいて上記暖房機の上記ONOFF制御を行う。また、上記暖房機の制御方法において、予め定められた直近の予測区間の上記複数の運転パターンについて上記温室の温度変化を予測するとよい。また、上記暖房機の制御方法において、上記ステップごとに又は複数の上記ステップごとに上記複数の運転パターンについて上記温度変化を予測し、上記最適運転パターンのうち上記温度変化の予測間隔分に基づいて上記ONOFF制御を行うとよい。また、上記暖房機の制御方法において、上記暖房機の制約条件を満たす複数の運転パターンについて上記温室の温度変化を予測するとよい。また、上記暖房機の制御方法において、使用頻度に応じて限定された複数の運転パターンについて上記温室の温度変化を予測するとよい。また、上記暖房機の制御方法において、上記暖房機は、木質燃料を燃焼させることで上記温室を暖房するようにしてもよい。また、上記暖房機の制御方法において、予測された上記温度変化の誤差と予め定められた燃費の重みとを含む評価関数に基づいて、上記温室の温度が目標温度に近づくように且つ燃費の悪化を抑えるように、上記複数の運転パターンの中から上記最適運転パターンを選択するようにしてもよい。本発明の暖房機は、温室に配置され、制御部によりONOFF制御される暖房機において、上記制御部は、制御タイミングとなる各ステップでのON又はOFFの組み合わせである複数の運転パターンについて、上記温室の物理モデルに基づいて上記温室の温度変化を予測し、予測された上記温度変化に基づいて、上記温室の温度が目標温度に近づくように、上記複数の運転パターンの中から最適運転パターンを選択し、上記最適運転パターンに基づいて上記暖房機の上記ONOFF制御を行う。本発明によれば、暖房機のONOFF制御を高精度に行うことができる。以下、本発明の実施の形態に係る、暖房機の制御方法及び暖房機について、図面を参照しながら説明する。<一実施の形態>図1は、本発明の一実施の形態に係る暖房機10を含む暖房システム1を示す概略構成図である。図1に示す暖房システム1は、暖房機10と、点火用灯油バーナ20と、押し込みファン30と、屋内サイロ40と、サイクロン集塵機50と、煙突60と、温室温度センサ91と、を温室100の内部に備える。また、暖房システム1は、灯油タンク70と、屋外貯留サイロ80と、外気温度センサ92と、を温室100の外部に備える。なお、温室温度センサ91は、便宜上図1の天井付近に二点鎖線で図示したが、例えば温室中央に配置して、計測温度が温室全体の温度の集中定数系として近似できるようにするとよい。また、外気温度センサ92も、便宜上図1の温室100の上方に二点鎖線で図示したが、例えば、平均的な外気温を検出できる位置に配置されればよい。暖房機10は、制御部11と、燃焼炉12と、温風吹出し口13とを有する。暖房機10は、制御部11によりONOFF制御され、木質ペレット201を燃焼炉12において燃焼させ、温風吹出し口13から温室100内に温風を供給する。木質ペレット201は、屋外貯留サイロ80から屋内サイロ40を介して暖房機10に供給される。また、木質ペレット201は、点火用灯油バーナ20により点火され、押し込みファン30により空気が暖房機10内に取り込まれた状態で燃焼する。木質ペレット201の点火は、灯油タンク70から供給される灯油が点火用灯油バーナ20で点火し、この点火により発生した火炎によって上記木質ペレット201を燃焼させることによって成される。なお、暖房機10から排出される排気ガス202は、サイクロン集塵機50によって燃焼灰203が取り除かれた後、誘引ファン204によって煙突60に送られ、温室100の外部に排出される。図2は、本実施の形態に係る暖房機の制御方法を説明するための温室温度モデル予測制御を示す概略フローチャートである。図2において、予測区間となるステップ数をP、後述する制約条件をC、目標温度をTref、誤差及び入力の重みをQ,R,制御時間をN、規定の繰り返し回数をDとする。なお、予測区間は、温室の物理モデルに基づいて温室の温度変化の予測をする直近の区間であり、例えば12ステップ1ステップ20秒などと予め定められる。予測区間の長さは適宜決定されればよいが、運転を行う全ての区間ではなく、その一部とすることで制御の複雑化を抑えることができる。まず、予測区間Pのそれぞれのステップ数においてON又はOFFの2値が与えられるため、2Pの入力列を生成する工程S1。そして、上記入力列の中から運転制約条件Cを満たす入力列を選定する工程S2。上記の予測区間Pについての入力列2Pの全ての運転パターン制御タイミングとなる各ステップでのON又はOFFの組み合わせに対して評価関数Jの値を求めると、その計算量は膨大なものとなる傾向がある。暖房機の制御アルゴリズムは安価な組込システム上で実現するのが望ましく、この要求からは出来るだけ計算量が小さい方が望ましい。ここで、暖房機に制約条件Cが発生する理由について簡単に述べる。図1に示す暖房システム1に運転命令が印加されると、屋内サイロ40から暖房機10の燃焼炉12に一定量の木質ペレット201が供給される。次いで、木質ペレット201に灯油タンク70から少量の灯油が散布されて点火用灯油バーナ20により着火が行われる。従って、暖房機10から定常的な熱供給が行われるようになるためには一定の時間が必要となるこれが木質ペレット201焚きの暖房機10が長時定数を有する理由でもある。よって、暖房機10のONOFF信号の切換えが頻繁に行われると、木質ペレット201の十分な燃焼が確保されないうちに消火が行われ、同時にまた種火を作るための灯油の消費量も増大し、燃費が大きく低下することになる。これを回避するためには、ON信号あるいはOFF信号を印加した後には、一定時間は運転信号を反転できない時間を設定するのが合理的であり、これが暖房機10の制約条件Cとなる。まず、予測を行う予測区間P、制約条件C、評価関数Jここでは予測区間Pにおける温度誤差の自乗和を与える。例として、予測区間P=3、C=2暖房機立上り/立下りとも同一とし、現時点での暖房機への入力が1=ONから0=OFFに初めて切換った場合を考える。考えられる全てのONOFFの入力列は23=8通り存在するが、図17に示す場合、1又は0が連続しない5パターンは採用できないため、運転パターンとして採用できるものは3通りしか存在しないことが分かる。図17の8つの運転パターンのうち、No.1は暖房機を全く運転しないもの全てOFF0、No.5は初めの2ステップをOFF0として最後のみON1にするもの、最後はNo.7で、初めの1ステップのみOFF0にして終わり2ステップ分はON1にするものである。詳しくは後述するが、こうして選ばれた3つの運転パターンに対し、暖房機への制御入力、温室温度センサ91により計測された温室温度、及び、外気温度センサ92による外気温度の情報を基に、予測区間Pステップ分に亘る温室温度の予測値を順次求めることができる。この結果に対し、評価関数J目標温度Trefとの誤差の自乗値の予測区間Pステップ間における総和の値が最小になる最適運転パターンを採用する。即ち、予測された温度変化に基づいて、温室の温度が目標温度Trefに近づくように、最適運転パターンを採用し工程S3,工程S4、初めの1ステップ分のONOFF信号を最終的に実際の運転信号として決定する工程S5。なお、最適運転パターンは、他の選択条件を基に選択されるものでもよい。以上の動作をステップ毎に繰り返して暖房機の制御を行う。具体的には、1つずつ増える制御時間N工程S6が規定の繰り返し回数Dに到達するまで行う工程S7。なお、2Pの入力列の生成工程S1は、繰り返し行う必要がないため、運転制約条件Cを満たす入力列の選定工程S2以降の処理が繰り返されることになる。なお、複数の運転パターンについて温室の温度変化の予測は1ステップごとでなく、複数ステップ例えば、数ステップごとに行うようにしてもよい。その場合、選択された運転パターンのうち温度変化の予想を行う間隔である予測間隔分のONOFFに従って制御を行うことになる。この場合、誤差が大きくなるおそれがあるものの、制御の複雑化を抑えることができる。ところで、制御対象のモデル予測器を内包したこのような予測に基づく制御は、動作環境の変化に柔軟に対応できること、制約条件C等を事前に考慮できることが特長である。ただし、ここでは、外気温が予測区間P中に変化しないと仮定している。さらに、同一予測区間Pにおいてほぼ同じ評価関数Jの値を達成する運転パターンが複数得られているならば、その中から出来るだけ暖房機を運転しない運転パターン或いは暖房機の着火の回数が少ない運転パターンを最適運転パターンとして選ぶことで燃費が最も良いものを選ぶことが可能となる。従って、温度誤差及び燃費の両方を考慮した評価関数Jを選ぶことができるのが、メリットの一つである。以下、温室モデルの構築、温室モデル予測制御、及び、最適解の分布について詳述する。温室モデルの構築温室の熱平衡モデルは、空間分布を考慮しない集中定数系として近似し、下記の数式1に従うものとする。ここで、T:温室温度、ht:熱貫流率、S:温室の表面積、V:温室の体積、ρ:空気の密度、Cp:空気の定圧比熱、Tex:外気温、Q:ヒータからの供給熱量、である。上記の数式1の左辺は、温室温度の時間変化率にそれ自身の変化のしにくさρVCpを掛けたもので、熱エネルギの時間変化量、すなわち温室に供給されるパワーを示している。一方、数式1の右辺は、温室内から逃げる熱エネルギ第一項目:温室内外の温度差により失われる熱エネルギで、その係数は被覆材であるビニール等の熱の通り易さ熱貫流率に面積を乗じることで得られると暖房機により温室内に供給される熱エネルギとからなる。温室に供給されたパワーと温度差により逃げた分と供給された分の差とがバランスすることを利用して得られたのが上記の数式1である。さらに暖房機からの供給熱エネルギは、最大供給可能な熱エネルギ量Qmaxで正規化する大きさを1に揃えることで、熱の発生状態を示す変数ξおよび効率ηとの積で下記の数式2に表現することができる。ただし、ξは0熱の発生が0の状態から1熱の発生が最大値となっている状態の間の任意の値をとるものとする。この熱の発生状態ξは、運転信号の印加に応じてさらに下記数式3の一次遅れ系に従うものと仮定する。ただし、K:ヒータの時定数、uONOFF:暖房機への制御入力ON=1,OFF=0である。この近似は暖房機の運転データより妥当性が認められる。以上の物理法則及び仮定から、温室モデルは二次系の連続時間システムとして下記の数式4及び数式5のようにモデル化される。これら数式4及び数式5を一定のステップで離散化することで、離散時間系における下記の数式6及び数式7の差分モデルが得られる。これら数式6及び数式7を単位時間遅延演算子q1ある信号xtに対し、q1xx、すなわちこの演算子が作用すると時刻が1だけ遅れるで書換えると、下記の数式8及び数式9が得られる。これよりξの動特性を消去し、外気温と暖房機の運転信号から温室温度までの動特性を求めると、伝達関数表現では以下の結果が得られる。この数式10の左式は、右式の伝達関数Gq1で表され、現時刻kにおける温室温度が、時刻k-1及びk-2における温室温度及び外気温、並びに、時刻k-2における暖房機の運転信号を基に決定されることを示している。従って、事前に得られている各物理定数を代入した後にシステムを離散化すれば、理論上は温室温度の数学モデルが得られることになるが、その結果を用いて求めた温室温度と実際の温度計測値とがかなり異なるのが現状である。この要因としては、熱貫流率が気象条件等により変化すること、地熱の影響を無視していること、システム自体を集中定数系として近似していること、等が挙げられ、これらの不確かさを含んだシステムは、物理的法則のみによる微分方程式で高精度には表現できないことが原因と考えられる。従って他の策としては、実際に計測された温室の内外気温及び暖房機の運転信号から各パラメータを数値的に推定することが挙げられる。もし温室とほぼ同じ動特性を示す数学モデルが入手できれば、それを設計に利用する公称値モデルとして採用できることになる。数式10で求められたように、温室温度の動特性伝達関数Gq1では、下記の数式11のように仮定できる。ここで、温度データについては、太陽の熱輻射などにより十分に温室温度が上昇する日中は除き、例えば19時から翌日8時までの間を制御対象とする。さらに、n4sid法を用い、精度の向上を図る。n4sid法とは、部分空間同定法とも呼ばれており、入出力データを基に拡大可観測行列と呼ばれる行列を計算し、その特異値を用いてシステムの次元を決定し、状態量、そして状態方程式を与えるシステム行列の同定を行う方法である。この方法により、温室のシステムの状態空間モデルを導出する。ただし、同種の同定アルゴリズムは多数存在し、計算量および精度により適切なものを選べばよい。システム同定に用いる温室データは、炉床式暖房機のものである。実際の温室データ20日分のうち、データとして信頼できると判断した温室データについてn4sid法によるシステム同定により温室パラメータを求め、公称値モデルを決定する。図3にシステム同定により得られた、伝達関数の係数のパラメータを示す。図3のPolesの列の値は、20日分それぞれの温室データから求めたシステムの極の絶対値である。これらの絶対値の値は、すべて1より小さくなっていることから、20日分のデータがすべてシステムとして安定であるといえる。よって、図3の20日分の温室データから平均をとり公称値モデルを決定することとする。ただし、パラメータのばらつきを勘案するために、その分散も考慮した。さらに、外気温を厳しい条件である3°Cにおいて温室温度の最終値が目標温度の20°Cに達するように各パラメータを調整した。求められたシステムの伝達関数を下記の数式12に示す。次に、実際の温室データの外気温と暖房機のONOFFの入力データとをn4sid法により求めた温室モデルに印加し、モデル出力の温度と実際の温室で測定された温度との比較を行う。図4に、n4sid法によるシステム同定により得られた温室モデルと公称値モデル温室温度の比較の例を示す。図4には、制御タイミングであるステップごとの温室温度greenhousetempを点線で、公称値モデルの出力nominaltempを実線で、外気温を破線で示した。図4では、温室の基礎計測データの中で、外気温が最も低く、温室温度及び公称値モデルの出力は、外気温よりも高い温度で、暖房機のONOFFに対応して上下に変動している。図4において、パラメータ調整後の公称値モデルと実際の温室温度との平均誤差率は2.72であった。一方、実際の計測データを基に同定を行ったモデルのパラメータが上記の数式1及び数式3でのそれに近いものであれば、同定によりパラメータを求める作業は不要となり、数式1及び数式3を基に公称値モデルを決定できることになる。以下では、物理モデルのパラメータがどの程度制御系の設計に利用可能かを考察する。温室の仕様から、放熱係数ht=6.7[W]、温室の表面積S=519.9[m2]、温室の体積V=678[m3]、空気の密度ρ=1.27[kgm3]、空気の定圧比熱Cp=1.0×103[J]であるため、物理モデルの上記数式1より、下記の数式13及び数式14が導き出される。一方、n4sid法によるシステム同定によって得られた公称値パラメータより、下記の数式15、数式16及び数式17が導き出される。この結果に基づく物理モデルと公称値モデルとの比較を図5に示す。まず、公称値パラメータから計算したα、γの値から求めた放熱係数htを、それぞれht,αとht,γとする。次にβの値から求めたヒータの効率をηβとし、物理モデルのパラメータと比較を行った。公称値モデルnominalvalueを基準とし物理モデルのパラメータphysicalvalueと比較すると、ht,α及びηβのパラメータについては物理モデルと公称値モデルとの対応ratioofparameterは約1.5~2倍程度になっていることがわかる。一方、ht,γのパラメータについては、外気温から温室温度への影響度合いを表すγが非常に小さな値となった。これは、同定計算においてはパラメータγが影響の小さな定値外乱として扱われてしまったことによるものであると考えられる。また、同定により得られた公称値モデルの傾向として暖房機の能力が低めに計算されていることが同定結果のデータ解析により明らかとなった。これは、温室の空間分布を考慮していないが、実際は暖房機からの熱エネルギーは空間分布を介して温室を温めるため、集中定数系を仮定した場合よりも同定結果は小さい値となる。つまり、見かけの暖房機能力が低めに計算されるため、外気温の影響を受けやすいシステムとして同定されたと考えられる。このことから、システム同定により得られた公称値モデルと理論から導いた物理モデルとの間には無視できないパラメータの誤差があることが分かった。実際に物理モデルのパラメータと温室データを用いて温室温度の再現性の確認を行ったが、実際の温室温度には全く追従しないという結果を得た。よって、物理モデルのパラメータをそのまま制御系の設計に用いることは難しいと考えられる。以降では、上記の数式12を用いる。温室モデル予測制御次に、温室モデル予測制御について説明する。ここでは、温室温度制御のシミュレーションを行った結果を述べる。シミュレーションは、1制御対象の温室とモデル予測制御系の予測器が一致している場合、及び、2α,β,γをシミュレーション中に強制的に変化させた場合について考察する。ただし、いずれも目標温度周り±1°Cの範囲での制御を目的とした例である。シミュレーションを行う際の予測器は、上述の公称値モデルを採用する。1制御対象の温室とモデル予測制御系の予測器が一致している場合まず、制御対象である温室と予測器のパラメータが一致している理想的な場合についてシミュレーションを行う。ただし、このシミュレーションでは、予測区間Pを12ステップ、運転制約条件CはON→OFF及びOFF→ONともに5ステップ、1ステップを20秒として計算を行った。4通りの異なる外気温条件におけるシミュレーション結果を図6~図9に示す。ここで、使用するこれら4通りの外気温は、20通りの温室データの中から夜間の外気温の変動幅が最も小さい場合、外気温の変動幅が最も大きい場合、外気温の平均値が最も高い場合、外気温の平均値が最も低い場合を選んだ。図6に示すように外気温の変動幅が最も小さい場合約7°C~約9°Cでほぼ安定のシミュレーション結果では、目標温度である20°Cに対する温室温度greenhousetempの平均誤差率は0.18で、-0.07~0.08°Cの範囲で温度制御可能であることがわかった。図7に示すように外気温の変動幅が最も大きい場合約7°C~約9.5°Cで上下動のシミュレーション結果では、目標温度である20°Cに対する温室温度greenhousetempの平均誤差率は0.15で、-0.06~0.06°Cの範囲で温度制御可能であることがわかった。これは、温度の高低差は大きいが急激な変動はなく比較的穏やかな温度変化であったことから予測器の予測出力の計算にほとんど影響しないためと考えられる。図8に示すように外気温の平均値が最も高い場合約11°C~約13.5°Cのシミュレーション結果では、外気温が高いため目標温度である20°Cに達するまでの時間は早いが、制御が難しく過渡反応が出ている。定常状態では、目標温度に対する温室温度greenhousetempの平均誤差率は0.21で、-0.08~0.08°Cの範囲で温度制御可能であることがわかった。図9に示すように外気温の平均値が最も低い場合約3.5°C~約6°Cのシミュレーション結果では、目標温度である20°Cに対する温室温度greenhousetempの平均誤差率は0.27で、-0.10~0.08°Cの範囲で温度制御可能であることがわかった。2α,β,γをシミュレーション中に変化させた場合実際の温室温度制御では、温室のシステムのパラメータは降雨や風等により、放熱係数が大きく変動する。このような状況を踏まえて、パラメータをシミュレーション中に強制的に変化させた場合について制御性を確認する。このパラメータ変化は実際の現象を考慮し、図10及び図11に示すように斜線部分の400Stepから600Stepにかけて徐々に変化を与えることで、設定したパラメータの値に近づけている。また、外気温は平均的な日のデータを用いている。このような変動のもとでも温度制御が高精度に行われるのであれば、実用性はさらに高くなる。まず、連続時間系でαパラメータに20の変化を与えたときの結果を図10に示す。図10のシミュレーションでは目標温度である20°Cに対する温室温度greenhousetempの温度誤差が0.35°Cで、平均誤差率1.39で温度制御が行えている。次に、連続時間系でγパラメータに20の変化を与えたときの結果を図11に示す。図11のシミュレーションでは目標温度である20°Cに対する温室温度greenhousetempの温度誤差が0.15°Cで、平均誤差率0.35で温度制御できている。図10及び図11に示す2通りのシミュレーション結果から、温室と予測器の公称値モデルの誤差を考慮することによりモデル化誤差の補償を行っている。これにより、モデル化誤差はある程度補償され、それにより目標温度との誤差は±1°C以内に制御できている。しかし、各パラメータが変化することにより無視できない定常誤差が残る結果となった。このときに印加した入力の傾向をみると、パラメータを変化させていない場合と同じ入力パターンであった。これは、予測器による予測出力の計算は誤差の自乗和が最小となるように正確に行えているが、その計算により決定した入力を、パラメータを変化させた温室に印加すると、図6~図9に示すシミュレーションと同様に、温室のパラメータの変化量に振られて定常値が狂ってしまうことが考えられる。つまり、パラメータ変動に対するある程度のロバスト性は認められるものの、その大きさはそれほど大きなものでないことがわかる。これは、モデル予測制御においては、予測器の精度に制御性能が大きく依存することによる。この原因として、誤差の補償の効果はそれほど大きなものではないことが挙げられる。このことから、温室温度を目標温度に近づけるためには、誤差の大きさを適時確認し、それが大きい場合には再度制御対象の温室のパラメータを同定する適応同定を行うことで改善されると考えられる。最適解の分布次に、最適解の分布について説明する。上述の例では、予測区間Pにおける2P個の運転パターンの中から暖房機の運転制約条件Cを満たすものだけを抽出するステップを踏んでいる。しかし、実際には数値シミュレーションを行ううちに、ほとんどの場合、制約条件Cを満たした後の候補の中でも選ばれる運転パターンがほぼ決まったものになることが明らかとなった。そのため、評価関数の値を計算する運転パターンを初めから絞り込む使用頻度に応じて限定することができることが判明した。これを示したのが後述する図12~図15のグラフである。これらのシミュレーションでは予測区間P=12と選んでいるため、全ての運転パターンは4096通り存在する。図12~図15では横軸に時間ステップを、縦軸に全ての運転パターンをとり、各ステップで何番目の運転パターンが最適運転パターンとして採用されたかを示している。図12~図15より、選ばれる最適運転パターンは、ほぼ4グループグループ1~4に分類されることが分かり、これ以外を計算する必要性は極めて低いことが分かる。これを基に最適化計算を行うことで、計算量は従来の6割程度まで抑制できることが明らかとなった。ただし、この結果は外気温及び制約条件に依存するため、それらに応じて適宜抽出する必要があるといえる。以下では、温室モデル予測制御における最適解の分布を調べ、その経験則をもとに演算量の低減化を図れるかについて検討した結果を述べる。演算量の低減化を定量的に示すためシミュレーションの演算時間を計測し、さらに最適解の分布対策の有無による制御性についても確認する。このとき分布対策とは、入力パターンに偏りがある場合に、予測出力の計算範囲を偏りのある範囲に限定し演算量の低減化を図るものである。ただし、外気温は平均的な変化を示す日のものを用いている。まず、温室と予測器が一致している場合の最適入力の分布を図12に示す。図12の縦軸は各Stepにおける最適解の入力番号を示している。このとき、実際の制御実績は100Step程で温室温度が定常状態になっており、図中には太線で示した。図12の結果より、定常状態になったとき最適入力の分布には偏りがあり、ある決まった入力パターンのみが選択されていることが分かる。次に夜間外気温の平均値が低い場合図13、夜間外気温の平均値が高い場合図14、及び、連続時間系で予測器のパラメータの値を変えたもの図15について、シミュレーションを行い、最適入力分布を調べた。まず、図13及び図14に示すように夜間外気温の違いによって最適解の分布に大きな変化は見られなかった。また、それぞれパラメータを変化させたα図15、β及びγの3通りについても最適解の分布を調べたところ、この場合も温室と予測器が同じ場合に採用される入力パターンと強い類似性があることが分かった。そのため、これまでのシミュレーション結果を基に最適入力の分布範囲を1から500番目、951から1050番目、1951から2150番目、3051から4095番目の4区間に絞り、これらの限定した入力候補に対して予測出力の計算を行い、どの程度演算量が低減化されているかを調べた。また、その時の制御結果を図16に示す。制御結果は入力候補を限定する前と後で全く変わっていない。これは、限定した入力パターンの区間だけに実際の制御入力が集中していたため制御性には影響が出なかったためと考えられる。また、1080Step分のシミュレーションにかかる計算時間は、入力候補を限定する前に約17秒かかっていたが、入力候補限定後には約10秒に短縮された。今回のシミュレーションでは予測区間Pを12Stepと定めているため、その時点で入力パターンは4096通りに決まるただし、厳密には全てOFFの入力パターンを計算から除いているため実際は4095通りである。しかし、経験則をもとに入力候補を限定したことにより、図16のシミュレーションでは入力パターンが1845通りになっているため、予測出力の計算に要する時間が短縮されたことが考えられる。この結果から、入力の候補数によって演算時間が短縮されるため、入力の候補を限定する範囲をさらに小さくすることで演算量はさらに削減されると考えられる。以上説明した本実施の形態では、制御タイミングとなる各ステップでのON又はOFFの組み合わせである複数の運転パターンについて、温室100の物理モデルに基づいて温室100の温度変化を予測し例えば、図2の工程S3、予測された温度変化に基づいて、温室の温度が目標温度Trefに近づくように、複数の運転パターンの中から最適運転パターンを選択し、例えば、図2の工程S4、選択された最適運転パターンに基づいて暖房機10のONOFF制御を行う。そのため、例えば、目標温度Tref例えば、20°Cとの誤差の少ない最適運転パターンを用いて、温室100の温度を目標温度Trefに近づけることができる。よって、本実施の形態によれば、暖房機10のONOFF制御を高精度に行うことができる。また、本実施の形態では、予め定められた直近の予測区間Pの複数の運転パターンについて、温室100の温度変化を予測することで、制御の複雑化を抑えることができる。また、本実施の形態では、ステップごとに又は複数ステップごとに複数の運転パターンについて温度変化を予測し、予測間隔分に基づいてONOFF制御を行うことで、ONOFF制御を非常に高精度に行うことができる。また、本実施の形態では、暖房機10の制約条件Cを満たす複数の運転パターンについて温室100の温度変化を予測することで、制御の複雑化を抑えることができる。また、本実施の形態では、全運転パターンのうち使用頻度に応じて限定された複数の運転パターンについて温室100の温度変化を予測することで、制御の複雑化を抑えることができる。また、本実施の形態では、暖房機100は、木質ペレット木質燃料201を燃焼させて温室100を暖房することで、重油燃料のような大きな燃料価格の変動ひいては商品価格の変動を抑えることができる。また、二酸化炭素排出量を抑えることができる。また、本実施の形態の制御方法により、外気温変化、暖房機の性質などに依存しない温度制御性を得ることもできる。ところで、わが国の植物栽培用の温室の栽培面積は、5万haに及んでおり、その内44に相当する2万2千haで暖房機が使用されている。したがって、本実施の形態の暖房機の制御方法及び暖房機をそれらに適用することが可能である。なお、本実施の形態では、木質ペレット201を燃料として用いる暖房機10を例に説明したが、ONOFF制御される暖房機であれば、他の木質燃料、オイル、或いは、これら以外の燃料を用いる暖房機を採用してもよい。また、暖房機10としては、温風により温室100を暖房するものではなく、温水などにより温室100を暖房する暖房機を採用してもよい。また、本実施の形態では、暖房機10の制御をONOFF制御のみについて説明したが、燃料である木質ペレット201の供給量や、押し込みファン30による空気の取込み量などの他の制御も行うものであってもよい。また、本実施の形態では、予測区間Pを12ステップと設定したり、制御タイミングとなる1ステップを20秒と設定したりするなど、制御の具体例を示すために一例としての数値を適宜設定しているが、このような数値は、温室の環境、温室に求められる温度制御の精度などによって適宜変更されるものである。<他の実施の形態>本実施の形態では、モデル予測制御の概念による温室の温度予測を基にした最適制御手法の有効性について、実際の温室における実験により確認を行った。今回の試験結果より、モデル予測制御アルゴリズムの導入により、実際の暖房機に対してもONOFF信号の切換えに対する制約の影響を考慮しながら、目標温度周りでの最大誤差をpp値で1.2°C以下に制御可能であることが確認された。試験1及び2は、ネポン株式会社厚木事業所東ハウスにおいて行った。試験1は、2011年4月27日22:05~23:45強風に行った。試験2は、2011年4月28日3:50~5:30雨強風に行った。試験1及び2では、目標温度Trefは25°C、1ステップは60秒、予測区間Pは720秒12ステップ、暖房機立上り/立下りの制約条件Cは300秒5ステップである。なお、予測制御で重要となる制御対象の予測器モデルのモデルパラメータは、2011年4月22日~26日に厚木事業所東ハウスで事前計測した温室データに対するシステム同定結果より、上述の一実施の形態で述べた手法で下記の数式18が得られた。図18は、温度制御結果1及び比較例1の温度シミュレーション結果を示すグラフである。図18に示すように、温度制御結果1Proposedmethodでは、温度誤差はpp値で目標値周り1.1°Cないし1.2°Cに制御されており、非常に高い精度で目標を達成できていることが読み取れる。なお、図18には暖房機のONOFF信号も合わせてプロットしている。得られた温度制御結果は図20及び図21にまとめられている。図20に示すように、温度制御結果1においては、最大誤差が0.7°C、pp値が1.2°C、平均誤差が0.27°C、目標温度reftempが上記のとおり25°Cである。また、図18には従来の運転アルゴリズム比較例1で暖房機の運転を行った場合についてシミュレーションを行った結果Conventionalmethodも合わせて示している。ここで、比較例1は、目標温度に対し1°Cを暖房機のOFF温度、目標温度に対し-1°CをON温度に設定し、目標温度との温度差のみから暖房機の運転信号を決定する方法のことで、このときのシミュレーションには上記の温室モデルの動特性を実際の温室の動特性として採用して計算を行っている。この場合には、図21に示すように、最大誤差が2.99°C、pp値が5.41°C、平均誤差が1.58°C、目標温度が上記のとおり25°Cである。これは、目標温度の±1°CにONOFFの閾値を設定しても、温室温度の時定数が非常に大きいために誤差が大きくなったものと考えられる。このことは、外気温変化がより大きい条件の下では運転制約のために制御性能がさらに悪化する可能性のあることを示しているが、温度制御結果1ではこの影響に対しても運転制約を考慮しながら最適な暖房機の運転が行われるため、従来の運転アルゴリズムである比較例1よりも制御性能の悪化が小さく抑えられると期待される。図19は、温度制御結果1及び比較例2の温度シミュレーション結果を示すグラフである。図19には、上述の温度制御結果1と制御結果のグラフは図18と同一である。、従来の運転アルゴリズム比較例2で暖房機の運転を行った場合についてシミュレーションを行った結果Conventionalmethodと、を合わせて示している。比較例2は、目標温度に対し閾値を±2°C実際に稼働しているペレット焚き暖房機の標準的な設定に変更した場合の方法である。この場合には、温度制御性能はさらに低下し、図22に示すように、最大誤差が3.75°C、pp値が7.16°C、平均誤差が2.05°C、目標温度が上記のとおり25°Cである。比較例2では、比較例1よりも閾値の絶対値が大きくなっているため、pp値などが比較例1の値よりも大きくなっている。図23は、温度制御結果2及び比較例1の温度シミュレーション結果を示すグラフである。温度制御結果2は、提案手法に燃費も考慮するアルゴリズムを追加した場合の温度制御実験の結果である。具体的には、下記の数式19で示されるkステップ目の評価関数Jk、誤差eの二乗和と暖房機のON回数の総和に重みζを掛けたものとを加えて予測区間長12ステップに亘って考慮したもので構成されている。これを最小化するような制御入力を毎ステップごとに選ぶことを目的とする。これにより、複数の運転パターンの中から最適運転パターンを選択する際に、温室の温度が目標温度に近づくように且つ燃費の悪化を抑えるように選択を行うことができる。ここで、ζは暖房機を燃費させている総和区間に対する重み燃費の重みである。試験2における各設定パラメータは試験1の場合と同一に設定し、事前の温度制御シミュレーション結果から燃費項の重みζを0.5に設定している。その結果、定常状態になったと考えられる100ステップ分を比較すると、試験1と試験2とでは、暖房機を焚いている合計ステップ数は50ステップと49ステップとであり大きな差異はないが、OFFからONへの切換回数が9回から7回に減少することが確認された。具体的には、図20に示すように、温度制御結果2においては、最大誤差が0.7°C、pp値が1.1°C、平均誤差が0.32°C、目標温度が上記のとおり25°Cである。この結果は、燃費項燃費の重みを考慮したことの効果と考えられる。従って、温度制御性能を許容する範囲内で緩めることにより、温度誤差および燃費との最適化を図れることが明らかとなった。なお、温度制御結果2の条件下における比較例1では、図21に示すように、最大誤差が2.78°C、pp値が5.27°C、平均誤差が1.57°C、目標温度が上記のとおり25°Cである。図24には、上述の温度制御結果2と制御結果のグラフは図23と同一である。、従来の運転アルゴリズム比較例2で暖房機の運転を行った場合についてシミュレーションを行った結果と、を合わせて示している。温度制御結果2の条件下における比較例2では、図22に示すように、最大誤差が3.65°C、pp値が7.10°C、平均誤差が2.00°C、目標温度が上記のとおり25°Cである。本実施の形態では、予測された温度変化の誤差eと予め定められた燃費の重みζとを含む評価関数Jに基づいて、温室の温度が目標温度に近づくように且つ燃費の悪化を抑えるように、複数の運転パターンの中から最適運転パターンが選択される。そのため、暖房機のONOFF制御を高精度に行うことと燃費を向上させることとを両立させることができる。上述の本実施の形態の温度制御試験では、事前に計測した温室データからシステム同定を行い、それにより得られた温室モデルを予測器に採用して実験を行った。さらに温度制御のみでなく、燃費評価を加えた場合についても検討を行った。なお、2日間を通して夜間外気温変動は小さかったものの、両日とも試験中に強い風雨があり、温室の放熱係数に影響を与えたことが予測される。得られた温度制御試験結果からは、サンプリングタイムを一実施の形態の20秒=現状の温度計測周期から60秒に延長した場合でも、モデル予測制御手法を応用することで暖房機の運転制約条件を考慮しながら、目標とする温度周り±1°Cに制御可能であることが確認された。また、暖房機の燃費評価も考慮することにより、暖房機の切換回数の抑制が可能であることも確認された。なお、モデル予測制御は予測器の精度が大きく制御精度に影響しており、今後予測器内に適応同定機構を導入することでさらに精度の高い温度制御性能が得られると期待される。さらに適応機構の導入により、ペレットの種類の違いホワイト、バーク等による最大熱量のばらつき等にも対応可能である。"