2変数入力のシステム同定法によるビニールハウス内の気温・飽差モデル―春秋期におけるモデルの短期間構築と評価―

(1)

原著論文

2 変数入力のシステム同定法によるビニールハウス内の気温・

飽差モデル―春秋期におけるモデルの短期間構築と評価―

中山信 *

1）

_，高田拓

1）

_{，木村竜士}

1）

_，岡宏一

2） 1）高知工業高等専門学校〒 783-8508 高知県南国市物部乙 200-1 2）高知工科大学システム工学群〒 782-8502 高知県香美市土佐山田町宮ノ口 185 要旨ビニールハウス内の気温や飽差のモデル構築は，ハウス環境を予測して制御する場合に重要である．本研究では，面積 100 m2_{程度の小型ハウスですぐに制御環境を整えられるように，入力変数の数を抑えた単純なシステム} 同定法による気温・飽差モデルの構築を行った．ピーマンを栽培しているビニールハウスにて，2018 年秋と 2019 年春の無加温の時期に計測したデータを利用し，代表的な線形 ARX（Auto-Regressive eXogenous）モデルと OE （Output Error）モデルによる気温・飽差モデルを作成し，評価した．構築モデルは，ハウス外気温とハウス内日射量の 2 変数入力からハウス内気温を出力する気温モデルと，ハウス外飽差とハウス内日射量の入力からハウス内飽差を出力する飽差モデルとした．評価の結果，モデル構築に必要なシステム同定期間は 5 日間程度で十分であり，出力気温（飽差）を実測値と比較すると，春期には 2 週間程度まで Fit 率 70（75）% 程度の高い精度を維持できるが，秋期には Fit 率 70（80）% 程度を維持できるのは 1 週間以内であった．本手法では，モデルの入力変数が少ないため，設置センサ数を少なくでき，センサ設置後 1 週間以内にはモデル構築を始められることが示唆された．キーワード

ビニールハウス，気温，飽差，システム同定，ARX（Auto-Regressive eXogenous）モデル，OE（Output Error）モデル

はじめに

国内では，ビニールハウスなどを利用する施設栽培の生産率を高めることが課題となっている．天候に左右されずに，野菜等の安定供給を確保するためには，環境制御装置を導入した温室の割合を高め，生産性を向上させることが重要である（農林水産省 2017）．今後，更なる生産性の向上やエネルギーコスト削減のために，ハウスのモデルベース制御，特に，モデル予測制御を行える安価な環境制御装置が求められる．ハウス内の環境変数の値を予測し，制御するモデル予測制御を実施するためには，ハウス内外の環境変化を多入力多出力で統合的に扱え，全ての制御装置の制約条件を考慮することが必要である．ただし，日本のハウス環境はその土地特有の気象条件に大きく依存しており，ハウスの大きさ，形状，素材，制御装置，栽培植物などが各ハウスで異なる．ハウス環境全てを網羅する物理モデルを構築するには，個々のハウス毎に多くの物理定数を取得する必要がある．しかしながら，日本国内に多く点在する小型のハウスを対象とした物理モデルを一つ一つ構築することは現実的ではない．近年，比較的安価なマイコンやセンサ類であっても，高精度なデータを取得できるようになり，農業分野においても，環境計測装置の製作・設置が進んできた（中野ら 2018，農研機構西日本農業研究センター 2020）．特に，植物の生長に直接関係する気温や湿度などは，環境計測で必ず取得されるデータであり，気温モデルや湿度モデルの構築が行われている（Rodríguez et al. 2015）．しかし，植物生 * Corresponding Author E-mail: [email protected]

(2)

長の観点から，一般的に計測される相対湿度ではなく，水蒸気が飽和するまでに必要な水蒸気量である飽差が近年では利用される．飽差と植物の生育との関係を調べた研究は多く，飽差の制御が潅水の節約とともに，果実収量の増大につながるなどの報告が多くなされている（Leyva et al. 2015，Zhang et al. 2015，Zhang et al. 2017）．

環境条件が多様な面積 100 m2_{程度の小型ハウスでは，} システム同定法によるモデル構築が利用しやすい．システム同定法では，入力と出力データからその数値モデルを推定するため，単純で実際的な数値モデルが得られる．代表的なブラックボックスモデルである線形 ARX（Auto-Re-gressive eXogenous）モデルでは，ハウス外気温，ハウス外日射量，雲量や風速などを入力として，再現性の高い気温モデルが構築されている（Frausto et al. 2003，Nakayama et al. 2018）．さらに発展型のモデルとしてニューラルネットワーク（Neural Network）や遺伝的アルゴリズム（Genetic Algorithm）を利用した気温モデル構築が行われ，精度の向上がなされている（Frausto and Pieters 2004，Patil et al. 2008， Guzmán-Cruz et al. 2013）．これらのシステム同定モデルを用いて，気温（Pawlowski et al. 2011）や相対湿度（Outanoute et al. 2016）のモデル予測制御が行われている．より良いモデル構築のためには，入力変数の数を増やし，より複雑で高度なモデルを利用し，システム同定期間を増やすことが一般的である．しかしながら，旧来の小型ハウスなどでは，環境センサなどの設置状態が十分でない場合が多い．一方で，そのようなハウスでは，センサを設置してモデル予測制御を始めることで栽培環境を最適に保つことができ，収穫物の収量増や品質向上などにつながることが期待できる．センサシステムの導入を容易にするために，入力変数の数やシステム同定期間を抑えた精度の良いモデル構築が求められている．本研究では，小型ハウスですぐに制御環境を整えられるように，入力変数の数を抑えた単純なシステム同定法による気温・飽差モデルの構築を行う．自然換気のみで気温管理をする無加温の春期と秋期に，ピーマン栽培ハウスで得られた環境データを利用した．システム同定法により，低次離散時間伝達関数モデルとして表し，ハウス内気温やハウス内飽差を精度よく再現できる条件を調査し，典型的な小型ハウスでの利用を検討した．

システム同定法

ビニールハウス内の気温と飽差のモデルとして，様々なモデルが利用されてきたが，代表的なシステム同定モデルである ARX モデルと OE（Output Error）モデルを用いる．出力である気温や飽差は時間に明示的に依存せず，入力によってのみ出力が決定されると考えられるため，線形時不変（LTI：Linear Time-Invariant）の系として扱う（足立 2009，室井ら 2015）．離散時間 LTI システム線形時不変モデルにおいて，連続時間系の出力 y と入力 u の関係は（1）式の常微分方程式で表わされる．an と bm は入力と出力のモデルパラメータであり，n は出力次数， m は入力次数である． dn_y(t) dtn +an−1 dn−1_y(t) dtn−1 +· · · + a1 dy(t) dt +a0y(t) =bmd m_u(t) dtm +bm−1 dm−1_u(t) dtm−1 +· · · + b1 du(t) dt +b0u(t) （1）また，離散時間系の出力 y と入力 u の関係は，（2）式で表される．G は離散時間伝達関数とし，k をサンプリング周期毎の時刻（整数）として，q−1_{y(k) = y(k − 1)}_となるようなシフトオペレータ q を導入した．出力 y(k) と入力 u(k) を用いて，離散時間伝達関数 G を表すことができる． y(k) =b1q−1+b2q−2+· · · + bmq−m 1 + a1q−1+· · · + anq−n u(k) = G(q)u(k) （2）システム同定モデル離散時間 LTI システムでの出力である（2）式に，雑音項を加えると，雑音を考慮した離散時間 LTI システムの一般的な出力となる． y(k) = G(q)u(k) + H(q)w(k) （3）ここで，H(q) は雑音モデルであり，w(k) は平均値 0 の白色雑音とする．（3）式がシステム同定モデルの式となる．（2）式の分子分母パラメータを θ =[a, b] と置き，以下の評価関数 J(θ) を最小にする θ を探索することで，G(q)，H(q) を同定できる． J(θ) = 1 N N k=1 [H−1_{(q){y(k) − G(q)u(k)}]}2 _（4）ここで，N はシステム同定の際の全データ数である．図 1 に，2 つの代表的なシステム同定モデルの概要を，解析で用いる2変数入力系として記述した．ARXモデルは，式の誤差として雑音 w が記述される式誤差モデルである．エネルギー収支などから得られる第一原理モデル（物理モデル）に近いため，実時間予測や制御に適している．一方， OE モデルは，出力直前に雑音 w を考慮するため，雑音が出力に混入する出力誤差モデルである．各入力の出力への影響に特化したモデルであるため，時間応答シミュレーションに適している．

(3)

データと分析手法

ピーマンを栽培しているビニールハウスでの計測データを用いて，ハウス内気温およびハウス内飽差モデルを構築した．以下では，計測データ，およびシステム同定モデルへの入力値換算に必要な指標などについて述べる．ビニールハウス環境高知県安芸郡芸西村にある高知工科大学の実験ビニールハウス施設で計測を実施した．ハウスは南北棟，形状は幅 5.5 m，奥行 16 m，高さ 3.4 m のカマボコ形，外張り（PO フィルム）と内張り（ビニール）の二重構造であり， 12 月∼ 3 月の加温期を除いては自然換気のみの気温管理にてピーマン栽培が行われている（図 2 上）．図 2 下に 10 月中のハウス内の状態を示す．栽培品種はミオギ，裁植密度は畝間 1.6 m 株間 70 cm として，1 畝に 10 株，3 畝 30 株が植えられている．根元には潅水チューブを這わせてあり，タイマーバルブが日中に 4 回ほど 30 分間ずつ開くように設定されている．今回は，2018 年 8 月 27 日の定植後，10 月 15 日から 12 月 5 日，および 2019 年 4 月 20 日から 5 月 16 日の無加温期を計測期間として実験を行った．なお，12 月 1 日∼ 2 日の夜間に限り，ハウス内気温が 14 C 設定を下回った際に加温した．内張りと外張りの天窓は，巻き上げ式天窓自動開閉装置によりハウス内気温 27 C の設定にて連動して開閉されており，この時期の晴天時の日中は側面の巻き上げ式手窓や入口が概ね全開，夜間は全閉となっていることが多かった．データロガーとして GRAPHTEC 社の GL100 を設置し，ハウス内外の温湿度センサ（GS-TH），ハウス内照度センサ（GS-LXUV）のデータを 1 分値で取得した．ハウス内外の温湿度センサは，それぞれ日射を避けるための通風筒に入れられ，AC ファンにより強制通風が行われた．また，ハウス内の温湿度センサはハウス中央 150 cm の高さに，ハウス外の温湿度センサはハウス入口前に 50 cm の高さに設置した．照度センサはハウス内の中央 180 cm の高さに設置して，ピーマン群落の影に掛からないようにした．照度センサをハウス内に設置した理由は，ビニール越しに植物が受ける日射を正確に測るためである．内張りのビニールは 2018 年 8 月に新しく張り替えられたが，外張りの PO フィルムは 2009 年以来張り替えられていないため，やや汚れており，ハウス内照度［lx］はハウス外照度［lx］の約 70% であった．飽差「はじめに」で述べた通り，植物の生育は，相対湿度よりも飽差に強い影響を受ける．飽差［g/m3_{］とは，飽和} 水蒸気量［g/m3_{］と絶対湿度［g/m}3_{］との差であり，空} 気中にどれくらい水蒸気が入る余地があるかを示す指標である．施設栽培における一般的な適正飽差値は 3–6 g/m3 程度であり，葉の気孔が開き，蒸散や CO2の吸収が促され，光合成による生長が維持される（中野ら 2018）．飽差は，気温と相対湿度より算出する．気温 T［ C］での飽和水蒸気圧 e［hPa］は，Tetens の近似式より次式で与えられる．図 1 2 変数入力における ARX モデルと OE モデル図 2 ビニールハウス環境

(4)

e(T) = 6.1078 × 10T+237.37.5T （5）飽和水蒸気量 c［g/m3_{］は，理想気体の状態方程式により} 以下のようになる． c(T) = 217 ×_{T + 237.15}e(T) （6）飽差 HD（Humidity Deﬁcit）［g/m3_{］は，飽和水蒸気量}_c_［g/ m3_{］と相対湿度}_R_{［%］より計算できる．} HD(T, R) = c(T) ×100 − R₁₀₀ （7）本研究では，温湿度センサで取得した気温と相対湿度の計測値から，飽差値を算出した．照度と日射量農業分野では，照度［lx］より日射量［W/m2_］で植物の生育環境を評価するのが一般的である．日射量は放射照度と同義であり，放射照度［W/m2_{］は，単位面積あたり} の放射束［W］である．放射束［W］は，単位時間当たりに放射する光のエネルギーに相当する．一方，照度［lx=lm/ m2_{］は，単位面積あたりの光束［lm］である．放射束［W］} に標準比視感度［lm/W］を乗じて可視光波長領域で積分することで光束［lm］が求まる．本研究では，照度と日射量（放射照度）の変換式として，昼間の自然光における波長域 400–700 nm の換算値を用いた（Thimijan and Heins 1983）． 1 lx = 0.00405 W/m2 _（8）モデル評価の指標システム同定モデルの評価のため，シミュレーションの出力値 ˆy と実測値 y との比較を行った．出力値と実測値の関係を見るために，相関係数（Correlation Coeﬃcient）を利用した．また，モデルの適合度である Fit 率は，以下の式を使用した． Fit率= ⎛ ⎜⎜⎜⎜⎜ ⎜⎜⎜⎜⎜ ⎝1 − ��N k=1{y(k) − ˆy(k)}2 ��N k=1{y(k) − ¯y}2 ⎞ ⎟⎟⎟⎟⎟ ⎟⎟⎟⎟⎟ ⎠· 100% （9）ここで ¯y は実測値 y の平均値である．（9）式の右辺第二項は， RMSE（Root Mean Squared Error）を測定値の標準偏差で規格化したものであり，NRMSE（Normalized RMSE）と呼ばれる．測定データ図 3 と図 4 に計測データを示す．センサで取得されるのは，ハウス内照度，ハウス内気温，ハウス外気温，ハウス内相対湿度，ハウス外相対湿度である．照度，気温，相対湿度は日変化を示しているが，1 日の中での最高値と最低値の値は，天候の状態により影響を受けている．照度の変動の仕方から，典型的な天候と照度変化の関係を確認できるが，日中の照度値が極端に小さい日は雨天時であり，日中の照度値が大きく変動している日は曇天時であり，日中の照度値が大きく変化が小さい日が晴天時である．気温や相対湿度の変動の仕方は，照度に大きく影響を受けている．全期間でハウス内気温はハウス外気温より，日中は 4–6 C，夜間は 1–2 C 高い傾向を示しており，秋期の方が夜間の気温差は大きい．春期には，ハウス内の相対湿度は比較的ハウス外湿度と同じ傾向であり，秋期の日中には，ハウス内相対湿度はハウス外相対湿度よりも比較的高い傾向を示している．秋期の照度や気温などは期間を通じて徐々に小さくなっており，季節が移り変わっていることを示している．モデル入力の際は，計測された気温と相対湿度から飽差を算出し，照度を日射量に変換して利用した．分析方法モデルとしては，気温モデルと飽差モデルの 2 つを独立に構築した．ハウス内気温を出力するモデルの入力変数は，ハウス内日射量とハウス外気温とし，ハウス内飽差を出力するモデルの入力変数は，ハウス内日射量とハウス外飽差とした．モデル評価には，過去の測定データを利用してモデルを構築し，以降の測定データで評価するのが一般的である．まず，モデル構築に最適な測定データの期間（システム同定期間）の長さを調べ，構築したモデルがどの程度の期間有効に機能するのかを確認した．春期に関しては，最初に 2019 年 4 月 20 日から 4 月 30 日の期間のデータをモデル構築に用い，5 月 1 日から 5 月 15 日の期間を評価データとして扱った．モデル構築に利用するシステム同定期間として，システム同定期間の日数を変えてデータを選び，各々の期間毎にモデルを構築した．構築したモデルを用いて，評価期間のデータに対して，モデル出力値と測定値との間で，1 日毎に相関係数，RMSE， Fit 率を算出した．秋期に関しては，最初に 2018 年 10 月 15 日から 10 月 31 日の期間のデータをモデル構築に用い，11 月 1 日から 12 月 4 日の期間のデータを評価データとした．構築に利用するシステム同定期間の選び方やモデル評価の指標は春期と同等とした．両期間ともシステム同定モデルとしては， ARX モデルと OE モデルを使用した．

結果と考察

ビニールハウスの環境データである気温と飽差に関して，システム同定モデルにより気温モデルと飽差モデルを

(5)

構築し，測定データに基づいてモデル出力値を評価した結果を示す．また，モデル構築に利用するシステム同定期間による違いなどを比較し，実際のビニールハウスで利用しやすいモデル設定について検討した．気温・飽差モデルの構築と評価モデル構築の例を表 1 に示す．春期のデータに対して，システム同定期間を 4 月 20 日から 4 月 30 日とした場合である．気温モデルと飽差モデルに対して，各々 ARX モデルと OE モデルを作成することで，表 1 に示すようなモデル式を求めることができる．モデル構築時の Fit 率は気温モデルで 70% 程度，飽差モデルで 85% 程度であり， RMSE は気温モデルで 1.7 C 以下，飽差モデルで 0.57 g/m3 であった．実際には，システム同定期間の日数を変えながら，各々の場合にモデル式を導出した．構築したモデルに対して，評価期間の入力値を用いてモデル出力値を算出した．ここでは，異なる天候の状態を含んだ 5 月 12 日から 5 月 14 日（晴天，曇天，雨天）を評価期間とし，モデル構築に利用したシステム同定期間は 4 月 26 日から 4 月 30 日の 5 日間とした．図 5 に，モデルの入力値（内日射量，外気温，外飽差）と出力値（内気温，内飽差）を示す．気温モデルと飽差モデルの共通の入力であるハウス内日射量は，1 日の平均的な大きさが天候によって大きく変わり，突発的に太陽が出るような曇天時では晴天時よりも日射量の変動が激しい．各々のモデル入力値であるハウス外気温とハウス外飽差は，共通の入力値であるハウス内日射量と関連した日変化を示しているが，細部の変化はかなり異なる．気温モデルの出力であるハウス内気温は，ARX モデルと OE モデルとも計測値の変動を追従しており，計測値とモデル出力値との差である残差は， ARX モデルで±6 C 以内，OE モデルで±3 C 以内となった．春期におけるモデル評価とシステム同定期間との関係春期において，モデル構築のためのシステム同定期間を変更した場合，モデル出力値や評価指標がどのように変動するかを調べた．表 2 では，気温モデルの場合に，システム同定期間を 1 日ずつ増やした場合の評価指標の例を示す．評価期間の 5 月 1 日は共通とした．ARX モデルでは 4 月 30 日のみを用いてモデル構築した場合に，OE モデルでは 4 月 28 日から 4 月 30 日の 3 日間を用いてモデル 0 25000 50000 75000 100000 1 / 5 / 9 1 0 2 0 0 : 0 5 / 5 / 9 1 0 2 0 0 : 0 0 3 / 4 / 9 1 0 2 0 0 : 0 5 2 / 4 / 9 1 0 2 0 0 : 0 0 2 / 4 / 9 1 0 2 00:00 2019/5/150:00 照度 [lx ] 0 10 20 30 1 / 5 / 9 1 0 2 0 0 : 0 5 / 5 / 9 1 0 2 0 0 : 0 0 3 / 4 / 9 1 0 2 0 0 : 0 5 2 / 4 / 9 1 0 2 0 0 : 0 0 2 / 4 / 9 1 0 2 00:00 2019/5/150:00 気温 [℃ ] ハウス内ハウス外 30 50 70 90 1 / 5 / 9 1 0 2 0 0 : 0 5 / 5 / 9 1 0 2 0 0 : 0 0 3 / 4 / 9 1 0 2 0 0 : 0 5 2 / 4 / 9 1 0 2 0 0 : 0 0 2 / 4 / 9 1 0 2 00:00 2019/5/150:00 相対湿度 [%] ハウス内ハウス外 0 0.25 0.5 0.75 1 2019/4/20 2019/4/25 2019/4/30 2019/5/5 2019/5/10 2019/5/15 窓開度図 3 春期における照度，気温，相対湿度，窓開度の経時変化 2019 年 4 月 20 日 6:00∼ 5 月 16 日 6:00までのハウス内照度，ハウス外気温・ハウス内気温，ハウス外相対湿度・ハウス内相対湿度，窓開度である．

(6)

表 1 モデルパラメータの例

出力モデルモデル式 Fit 率［%］ RMSE

ハウス内気温 ARX

A(z)y(t) = B(z)u(t) + e(t)

A(z) = 1 − 1.98z−1₊_0.983z−2

B1(z) = 0.0220 − 0.0219z−1 B2(z) = 0.000394 − 0.000394z−1

67.5 1.71 C

OE

y(t) = B(z)_F(z)u(t) + e(t)

B1(z) = 0.950 B2(z) = 0.00276 F1(z) = 1 − 0.0989z−1 F2(z) = 1 − 0.879z−1 74.5 1.34 C ハウス内飽差 ARX

A(z)y(t) = B(z)u(t) + e(t)

A(z) = 1 − 1.535z−1₊_0.544z−2 B1(z) = 0.363 − 0.356z−1 B2(z) = 0.00108 − 0.000869z−1 84.6 0.571 g/m3 OE y(t) = B(z) F(z)u(t) + e(t) B1(z) = 0.366 B2(z) = 0.000930 F1(z) = 1 − 0.498z−1 F2(z) = 1 − 0.957z−1 85.1 0.552 g/m3 ※システム同定期間：4/20 ∼ 4/30，ハウス内気温モデル入力値：ハウス外気温 u1［ C］，ハウス内日射量 u2［W/m2］，ハウス内飽差モデル入力値：ハウス外飽差 u1［g/m3］，ハウス内日射量 u2［W/m2］ 0 25000 50000 75000 100000 8 1 0 2 0 0 : 0 4 1 / 1 1 / 8 1 0 2 0 0 : 0 4 / 1 1 / 8 1 0 2 0 0 : 0 5 2 / 0 1 / 8 1 0 2 0 0 : 0 5 1 / 0 1 / 8 1 0 2 /11/240:00 2018/12/40:00 照度 [lx ] 0 10 20 30 8 1 0 2 0 0 : 0 4 1 / 1 1 / 8 1 0 2 0 0 : 0 4 / 1 1 / 8 1 0 2 0 0 : 0 5 2 / 0 1 / 8 1 0 2 0 0 : 0 5 1 / 0 1 / 8 1 0 2 /11/240:00 2018/12/40:00 気温 [℃ ] ハウス内ハウス外 30 50 70 90 8 1 0 2 0 0 : 0 4 1 / 1 1 / 8 1 0 2 0 0 : 0 4 / 1 1 / 8 1 0 2 0 0 : 0 5 2 / 0 1 / 8 1 0 2 0 0 : 0 5 1 / 0 1 / 8 1 0 2 /11/240:00 2018/12/40:00 相対湿度 [% ] ハウス内ハウス外 0 0.25 0.5 0.75 1 2018/10/15 2018/10/25 2018/11/4 2018/11/14 2018/11/24 2018/12/4 窓開度図 4 秋期における照度，気温，相対湿度，窓開度の経時変化 2018年10月15日6:00∼12月5日18:00までのハウス内照度，ハウス外気温・ハウス内気温，ハウス外相対湿度・ハウス内相対湿度，窓開度である．

(7)

構築した場合に，最も Fit 率が高く，RMSE が小さく，相関係数が大きくなっているが，OE モデルの方が ARX モデルよりも Fit 率が高く，RMSE が小さい．図 3 に示しているが，照度，気温，相対湿度などの変動の仕方は 5 月 1 日と 4 月 30 日は比較的似ているが，4 月 29 日はやや異なり，4 月 28 日より以前は大きく異なっている．OE モデルは ARX モデルと比較して，同定期間でのモデル入出力値の変動を上手くモデル構築に反映できている．また，システム同定期間を増やすと，それぞれのモデルの Fit 率は下がり，RMSE は大きくなり，相関係数は小さくなった．システム同定期間のモデル入出力値の時系列変化が，評価期間の時系列変化の傾向と異なる日が多いほど，モデルは平均化された出力となる．表 2 のように評価期間を 5 月 1 日に定めた場合は，システム同定期間の日数が 1 日から 3 日 0 100 200 300 400 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 内日射量 [W /m 2] 10 15 20 25 30 35 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 外気温 [℃ ] 10 15 20 25 30 35 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 内気温 [℃ ] 計測値 ARXモデル OEモデル -6 -3 0 3 6 残差 [℃ ] ARXモデル OEモデル 0 5 10 15 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 外飽差 [g /m 3] 0 5 10 15 20 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 内飽差 [g /m 3] 計測値 ARXモデル OEモデル -4 -2 0 2 5/12 6:00 5/13 6:00 5/14 6:00 5/15 6:00 残差 [g /m 3] ARXモデル OEモデル図 5 気温・飽差モデル評価の例

(8)

で平均化したモデルが適切であることが分かる．ここで，表 2 で得られた気温モデルと飽差モデルの両方に関して，評価期間を 5 月 1 日から 1 日ずつ変更しながらモデル評価を行った．図 6 にはシステム同定期間の日数が 5 つの場合のモデルのみを示しているが，実際には 1 日から 11 日間の全ての場合について計算して調べた．システム同定期間の日数が 1 日や 3 日と少ない場合は，Fit 率などの変動が大きい．システム同定期間の日数が 5 日以上であれば，気温モデルで 75% 程度の Fit 率，飽差モデルで 80% 程度の Fit 率と非常によく計測値をモデル出力できている．ARX モデルと OE モデルでは，システム同定期間の日数が少ない場合に差があるが，5 日以上の場合は同等の Fit 率を示した．次に，システム同定期間を変更した場合のモデル出力値の変化を調べた．システム同定期間の日数は 5 日間に固定した．結果を図 7 に示すが，ARX モデルと OE モデルで大きな差はなかった．気温モデルではシステム同定期間が異なると同じ評価期間のデータに対しても Fit 率がややばらつくが，飽差モデルではシステム同定期間の変更によるばらつきは小さく，ばらつく場合も直近のデータを利用した場合の方が Fit 率は高くなった．一方で，4 月 28 日から 4 月 30 日の区間では，気温モデルと飽差モデルは Fit 率が極端に小さくなった．図 3 から推測できるが，4 月 28 日は曇天，4 月 29 日と 4 月 30 日は雨天であり，照度と気温，相対湿度の 1 日の変動の仕方が，他の期間と大きく異なっていた．そのため，この区間を評価した場合に大きく Fit 率が下がっている．今回利用している単純なモデルでは，モデル構築時に，モデル評価する日の変動パターンと似た日が含まれないと Fit 率が大きく下がる．今回の計測において，4 月 28 日から 4 月 30 日の期間は，測定を開始してから初めての雨天時であるため，最初のモデル構築時に，晴天と雨天が両方含まれるようにすれば，Fit 率は高くなると考えられる．結果として，春期の気温モデルと飽差モデルにおいては，概ね，良い Fit 率となっており，今回の単純なシステム同定モデルにより計測値に近いモデル出力が十分可能である．なお，図に示していないが，春期におけるシステム同定期間が 5 日間固定の場合，モデル出力値と計測値の RMSE は，気温モデルでは 0.70–3.0 C 程度で，飽差モデルでは 0.29–0.96 g/m3_{程度であった．} 秋期におけるモデル評価とシステム同定期間との関係秋期においても，同様にモデル出力値の評価指標について調べた．図 8 は，モデル構築のためのシステム同定期間の日数を変化させてモデル評価した例である．図にはシステム同定期間の日数が 6 つの場合のみ示しているが，実際には 1 日から 17 日間の全ての場合について計算して調べた．秋期においては，ARX モデルの方が全体的に Fit 率は低く，ばらつきが大きいが，特徴的な変化は OE モデルと似ている．以下では，OE モデルの変化の仕方を例に説明表 2 気温モデルの評価指標の例システム同定期間評価期間 Fit 率［%］ RMSE［ C］相関係数

期間日数 ARX OE ARX OE ARX OE

4/30 1 5/1 83.5 84.0 0.504 0.486 0.996 0.997 4/29–4/30 2 5/1 82.0 84.7 0.571 0.463 0.996 0.997 4/28–4/30 3 5/1 82.4 92.0 0.951 0.266 0.964 0.996 4/27–4/30 4 5/1 80.0 82.9 0.823 0.534 0.975 0.988 4/26–4/30 5 5/1 76.4 75.0 0.984 0.760 0.968 0.982 4/25–4/30 6 5/1 75.9 75.0 0.927 0.756 0.973 0.982 4/24–4/30 7 5/1 77.9 73.9 0.844 0.790 0.973 0.983 4/23–4/30 8 5/1 77.6 73.7 0.950 0.795 0.963 0.983 4/22–4/30 9 5/1 76.9 74.5 0.983 0.772 0.961 0.983 4/21–4/30 10 5/1 76.2 74.6 1.07 0.768 0.954 0.982 4/20–4/30 11 5/1 76.7 75.5 1.27 0.744 0.936 0.982 図 6 気温・飽差モデルでの Fit 率とシステム同定日数との関係（春期）

(9)

する．システム同定期間が 1 日のように少ない場合は，4 日以上の場合と比べて，Fit 率がばらつく．気温モデルでは，たまたま 1 日の場合の方が平均的に Fit 率は高いが，システム同定に選択した期間によっては，非常に悪くなることもあるため，モデル構築には 5 日間程度以上のデータ利用が望ましいと考えられる．評価期間全体での変動を見ると，気温モデルと飽差モデルは期間の後半になるにつれて，Fit 率が下がっている．図 4 で示したように，この期間が季節進行に対応しており，日中の照度や気温が徐々に下がっていることと関係している．季節進行のある時期に関しては，構築したモデルの有効期間は 1 週間程度と短いと考えられる．次に，システム同定期間を変更した場合のモデル出力値の評価指標を調べた．システム同定期間の日数は 5 日間に固定した．結果を図 9 に示すが，気温モデルと飽差モデルの両方で，評価期間の後半 11 月 14 日以降に Fit 率がばらついている．これは，システム同定期間が評価期間に近いデータでは，Fit 率が高くなっていることによる．秋期においては，Fit 率がある程度高い状態を維持できるのは，1 週間程度である．一方で，飽差モデルで顕著であるが，時折非常に Fit 率が悪くなる日がある．これらの日は，図 4 で日中の照度が小さい雨天時や曇天時に対応する．晴天時のみでモデル構築した場合に，雨天時を評価すると非常に Fit 率が下がる．一方で，雨天時や曇天時の変動の仕方は，日によって大きく変わる上，季節進行とともに変化しているようで，日にちが離れると変動の傾向を見いだせない．そのため，秋期に季節進行している際には，モデル構築に晴天時と雨天時を含んだとしても，日数が離れた期間の評価では Fit 率は高くならない．図 7 気温・飽差モデルでの Fit 率とシステム同定期間との関係（春期）図 8 気温・飽差モデルでの Fit 率とシステム同定日数との関係（秋期）

(10)

なお，図に示していないが，秋期におけるシステム同定期間が 5 日間固定の場合，モデル出力値と計測値の RMSE は，気温モデルでは 0.67–2.7 C 程度で，飽差モデルでは 0.17–1.2 g/m3_{程度であった．} 考察過去の ARX モデルなどを利用したシステム同定法による研究では，入力変数の数が 4 つ以上である場合が多く，天候や雲量などの状況を含んでいた（Frausto et al. 2003， Frausto and Perters 2004，Patil et al. 2008，Guzmán-Cruz et al. 2013）．本研究では，気温モデルと飽差モデルを独立させた上で，各モデルの入力変数を2つずつ（実質的には3変数）とかなり簡略化させている．それにも関わらず，特に春期においては十分高い Fit 率を出せており，小型ハウスなどの環境計測の設備が十分でないハウスであっても，容易にモデル構築を始められると考えている．過去の研究ではハウス外日射量を用いることが多いが，本研究ではハウス内日射量を使っており，天候や雲量などの情報以外に，ハウス被覆資材の情報も暗に含んだ状況となっている可能性がある．秋期においては，季節進行しながら晴天と雨天の周期的な入れ替わりがあるため，現状のモデルでは Fit 率が低い場合がある．例えば，1 年前の同じ時期の計測データから，晴天時と雨天時の日を選び，さらに直近数日の計測データを組み合わせて，モデル構築に利用するなどの方法を取ることで，秋期においても Fit 率を上げられる可能性がある．今後検討を行う予定である．

まとめ

本研究では，小型ビニールハウスですぐに制御環境を整えられるように，入力変数の数を抑えた単純なシステム同定法による気温・飽差モデルの構築を行った．ピーマン栽培のビニールハウスにて，2018 年秋と 2019 年春に計測したデータを利用し，代表的な線形 ARX モデルおよび OE モデルによる気温・飽差モデルを構築し，評価した．構築モデルは，ハウス外気温とハウス内日射量の 2 変数入力からハウス内気温を出力する気温モデルと，ハウス外飽差とハウス内日射量の入力からハウス内飽差を出力する飽差モデルとした．モデル構築に利用するデータは 5 日間程度で十分であり，出力気温（飽差）を実測値と比較すると，春期には 2 週間程度まで Fit 率 70（75）% 程度の高い精度を維持でき，秋期には 1 週間程度までは Fit 率 70（80）% 程度であるが季節の進行と共に Fit 率が低下した．気温モデルでは，モデル出力値と計測値の RMSE が春期で 0.70–3.0 C 程度，秋期で 0.67–2.7 C 程度であった．飽差モデルでは，実測値との RMSE が春期で 0.29–0.96 g/m3_{程度，秋期で 0.17–1.2} g/m3_{程度であった．一方で，秋期においてはモデル構築} に利用する計測データの組み合わせを考慮することで，Fit 率を上げられる可能性があり，今後検討を行う．結果として，ハウス内の気温・飽差のモデル構築において，センサ数を少なくし，センサ設置後 1 週間以内にはモデル構築を始めることができる．環境センサの設置状況が良くない小型ビニールハウスなどであっても，必要最低限のセンサシステムの導入はそれほど難しくないと考えられるため，利用しやすいことが示唆された．今後は，制御入力である窓開度やヒータ温度を含めたシステム同定モデルを作成し，夏期や冬期を含めた全期間で利用できるようにする．最終的には，ハウス環境のモデル予測制御システムの完成を目指す．図 9 気温・飽差モデルでの Fit 率とシステム同定期間との関係（秋期）

(11)

謝辞

ビニールハウス施設の利用に関しましては，高知県立高知工業高等学校北村晋助校長と，お母様の北村康子氏に大変お世話になりました．ここで深く御礼を申し上げます．本研究は JSPS 科研費 20K06332 の助成を受けたものです．引用文献足立修一（2009）システム同定の基礎，東京電機大学出版局，東京， 57–169．

Frausto, H. U., J. G. Pieters and J. M. Deltour (2003) Modelling greenhouse temperature by means of Auto Regressive Models, Biosystems Engi-neering, 84 (2): 147–157.

Frausto, H. U. and J. G. Pieters (2004) Modeling greenhouse temperature using system identiﬁcation by means of neural networks, Neurocom-puting, 56 (1): 423–428.

Guzmán-Cruz, R., E. Olvera-González, I. L. López-Cruz and R. Mon-toya-Zamora (2013) Modeling greenhouse air temperature using evolutionary algorithms in auto regressive models, African Journal of Agricultural Research, 8 (3): 251–261.

Leyva, R., C. Constán-Aguilar, E. Sánchez-Rodríguez, M. Romero-Gámez and T. Soriano (2015) Cooling systems in screenhouses: Eﬀect on mi-croclimate, productivity and plant response in a tomato crop, Biosyste-ms Engineering, 129: 100–111. 室井秀夫，竹下侑，足立修一（2015）時間領域におけるシステム同定結果の評価指標について，計測自動制御学会論文集，51 （10）：736–743．中野明正，安東赫，栗原弘樹（2018）ICT 農業の環境制御システム製作，誠文堂新光社，東京，24–43 & 100–105．

Nakayama, S., K. Oka and A. Harada (2018) Low cost modeling of green-house temperature and humidity deficit based on identification of transfer functions, SICE Annual Conference 2018, 1269–1272. 農研機構西日本農業研究センター（2020）簡易施設向け ICT 環境計測システム製作マニュアル，〈http://www.naro.aﬀrc.go.jp/publi city_report/publication/files/warc_man_ictpanfu20200330l.pdf〉， 2020 年 7 月 6 日参照．農林水産省（2017）施設園芸をめぐる情勢，〈https://www.maﬀ.go.jp /j/seisan/ryutu/engei/sisetsu/pdf/jyousei_1.pdf〉，2020 年 7 月 6 日参照．

Outanoute, M., A. Lachhab, A. Ed-dahhak, M. Guerbaoui, A. Selmani and B. Bouchikhi (2016) Synthesis of an optimal dynamic regulator based on Linear Quadratic Gaussian (LQG) for the control of the relative humidity under experimental greenhouse, International Journal of Electrical and Computer Engineering (IJECE), 6 (5): 2262–2273. Patil, S. L., H. J. Tantau and V. M. Salokhe (2008) Modelling of tropical

greenhouse temperature by auto regressive and neural network mod-els, Biosystems Engineering, 99: 423–431.

Pawlowski, A., J. L. Guzmán, F. Rodríguez, M. Berenguel and J. E. Normey-Rico (2011) Predictive control with disturbance forecasting for greenhouse diurnal temperature control, IFAC Proceedings, 44 (1): 1779–1784.

Rodríguez, F., M. Berenguel, J. L. Guzmán and A. Ramírez-Arias (2015) Modeling and control of greenhouse crop growth, Springer Interna-tional Publishing, Switzerland, 9–97.

Thimijan, R. W. and R. D. Heins (1983) Photometric, radiometric, and quantum light units of measure: a review of procedures for intercon-version, Horticultural Science, 18 (6): 818–822.

Zhang, D., Z. Zhang, J. Li, Y. Chang, Q. Du and T. Pan (2015) Regulation of vapor pressure deﬁcit by greenhouse micro-fog systems improved growth and productivity of tomato via enhancing photosynthesis during summer season, PLoS One, 10: e0133919.

Zhang, D., Q. Du, Z. Zhang, X. Jiao, X. Song and J. Li (2017) Vapour pres-sure deﬁcit control in relation to water transport and water productivity in greenhouse tomato production during summer, Scientiﬁc Reports, 7: 43461.

受付日 2020 年 8 月 18 日受理日 2020 年 11 月 29 日

(12)

Agricultural Information Research 30 (1) 2021. 1–12

Temperature and Humidity Deﬁcit Model for Greenhouses

Based on System Identiﬁcation Method with Two-Variable Input:

Quick Setup and Model Evaluation during Spring and Fall

Shin Nakayama*

1)

_{, Taku Takada}

1)

_{, Ryushi Kimura}

1)

_{and Koichi Oka}

2)

1) National Institute of Technology, Kochi College, 200-1 Monobe-Otsu, Nankoku City, 783-8508 Japan 2) Kochi University of Technology, 185 Miyanokuchi, Tosayamada, Kami City, 782-8502 Japan

Abstract

Building a model of temperature and humidity deficit (HD) is important for the predictive control of greenhouses. In this study, we developed temperature and HD models by using the system identification method with a small number of parameters for the quick setup of a control environment in a small greenhouse of about 100 m2_{. The temperature and HD models were created from} linear auto-regressive exogenous (ARX) and output error (OE) models. The modeling data were measured in a greenhouse where green peppers were grown unheated in the fall of 2018 and the spring of 2019. The greenhouse air temperature is output from the outside air temperature and the interior solar radiation. The greenhouse HD is output from the outside HD and the interior solar radiation. Data of 5 days are sufficient for building the model. Comparing the output temperature (and HD) with the measured val-ues, we found that the goodness of fit was highly accurate and could be maintained at 70% (80%) for 1 week in the fall and 70% (75%) for up to 2 weeks in the spring. Because the device needs few sensors and the model can be used within a week of installa-tion, the temperature and HD models can be used in small greenhouses without an extensive sensor network.

Keywords

greenhouse, temperature, humidity deficit, system identification, ARX (auto-regressive exogenous) model, OE (output error) model

* Corresponding Author E-mail: [email protected]