"外気冷房を含む複数の熱源機器を備えた空調熱源システムの運転を制御する熱源制御装置であって、所定の計測周期で計測された前記空調熱源システムの空調負荷および外気条件の計測値を時系列データとして格納する計測値データベースと、前記計測値の時系列データに基づいて、前記計測周期の所定回数先までの所定期間における前記空調負荷および前記外気条件の予測値を時系列データとして求める時系列予測部と、前記予測値の時系列データに基づいて、前記複数の熱源機器のモデルパラメータを推定するパラメータ推定部と、前記モデルパラメータに基づいて、所定の制御周期ごとの前記複数の熱源機器の運転状態を状態変数とし、少なくとも前記空調熱源システムの前記所定期間における運用コストの最小化を目的関数として定式化された混合整数計画問題を解き、少なくとも次の制御周期における前記状態変数を求める演算部と、前記演算部の演算結果に基づいて前記複数の熱源機器の運転を制御する制御部と、を有することを特徴とする熱源制御装置。前記演算部は、直前の制御周期における演算結果を利用して前記混合整数計画問題を解くことを特徴とする請求項1に記載の熱源制御装置。前記状態変数は、前記複数の熱源機器のそれぞれの停止または起動を表す0-1変数、および前記複数の熱源機器のそれぞれの出力を表す連続変数を含み、前記演算部は、分枝限定法を用いて前記混合整数計画問題を解き、前記分枝限定法においては、前記演算部は、前記0-1変数を0または1に固定する優先順位を決定するための指標値を求め、当該求めた指標値を、前記混合整数計画問題を解くのに利用するとともに、次の制御周期における前記指標値の初期値として利用することを特徴とする請求項2に記載の熱源制御装置。前記外気条件は、外気温度および外気湿度を含み、前記パラメータ推定部は、前記モデルパラメータが前記外気温度および前記外気湿度と対応付けられたパラメータテーブルを参照して前記モデルパラメータを求めることを特徴とする請求項1ないし請求項3の何れかに記載の熱源制御装置。請求項4に記載の熱源制御装置と、前記空調熱源システムと、前記外気温度および前記外気湿度を計測する温湿度計と、を備えることを特徴とする空調システム。前記複数の熱源機器は、冷凍機および冷却塔をさらにを含むことを特徴とする請求項5に記載の空調システム。外気冷房を含む複数の熱源機器を備えた空調熱源システムの運転を制御する制御部と、所定の計測周期で計測された前記空調熱源システムの空調負荷および外気条件の計測値を時系列データとして格納する記憶部と、を備えるコンピュータに、前記計測値の時系列データに基づいて、前記計測周期の所定回数先までの所定期間における前記空調負荷および前記外気条件の予測値を時系列データとして求める時系列予測処理と、前記予測値の時系列データに基づいて、前記複数の熱源機器のモデルパラメータを推定するパラメータ推定処理と、前記モデルパラメータに基づいて、所定の制御周期ごとの前記複数の熱源機器の運転状態を状態変数とし、少なくとも前記空調熱源システムの前記所定期間における運用コストの最小化を目的関数として定式化された混合整数計画問題を解き、少なくとも次の制御周期における前記状態変数を求め、前記制御部に入力する演算処理と、を実行させることを特徴とする熱源制御プログラム。前記外気条件は、外気温度および外気湿度を含み、前記記憶部は、前記モデルパラメータが前記外気温度および前記外気湿度と対応付けられたパラメータテーブルを予め格納しており、前記パラメータ推定処理は、前記パラメータテーブルを参照して前記モデルパラメータを求めることを特徴とする請求項7に記載の熱源制御プログラム。外気冷房を含む複数の熱源機器を備えた空調熱源システムの運転を制御する熱源制御方法であって、所定の計測周期で計測された前記空調熱源システムの空調負荷および外気条件の計測値の時系列データに基づいて、前記計測周期の所定回数先までの所定期間における前記空調負荷および前記外気条件の予測値を時系列データとして求め、前記予測値の時系列データに基づいて、前記複数の熱源機器のモデルパラメータを推定し、前記モデルパラメータに基づいて、所定の制御周期ごとの前記複数の熱源機器の運転状態を状態変数とし、少なくとも前記空調熱源システムの前記所定期間における運用コストの最小化を目的関数として定式化された混合整数計画問題を解き、少なくとも次の制御周期における前記状態変数を求め、求めた前記状態変数に基づいて前記複数の熱源機器の運転を制御することを特徴とする熱源制御方法。近年、さらなる省エネルギー化を図るため、熱源機器として、自然エネルギーを活用する外気冷房を備えた空調熱源システムも知られている。しかしながら、このような空調熱源システムでは、外気冷房の冷却能力は、温度や湿度などの外気条件から直接影響を受け、大きく変動してしまう。さらに、外気条件の変動に追従するため、熱源機器の運転制御を行う制御周期をより短くする必要がある。そのため、特許文献1の熱源運転支援制御方法を、外気冷房を備えた空調熱源システムに適用することは困難である。前述した課題を解決する主たる本発明は、外気冷房を含む複数の熱源機器を備えた空調熱源システムの運転を制御する熱源制御装置であって、所定の計測周期で計測された前記空調熱源システムの空調負荷および外気条件の計測値を時系列データとして格納する計測値データベースと、前記計測値の時系列データに基づいて、前記計測周期の所定回数先までの所定期間における前記空調負荷および前記外気条件の予測値を時系列データとして求める時系列予測部と、前記予測値の時系列データに基づいて、前記複数の熱源機器のモデルパラメータを推定するパラメータ推定部と、前記モデルパラメータに基づいて、所定の制御周期ごとの前記複数の熱源機器の運転状態を状態変数とし、少なくとも前記空調熱源システムの前記所定期間における運用コストの最小化を目的関数として定式化された混合整数計画問題を解き、少なくとも次の制御周期における前記状態変数を求める演算部と、前記演算部の演算結果に基づいて前記複数の熱源機器の運転を制御する制御部と、を有することを特徴とする熱源制御装置である。本発明の他の特徴については、添付図面及び本明細書の記載により明らかとなる。本発明によれば、外気条件を考慮して混合整数計画法を用いることにより、外気冷房を備えた空調熱源システムの運転制御を的確に行うことができる。本明細書および添付図面の記載により、少なくとも以下の事項が明らかとなる。===空調システム全体の構成の概略および熱源制御装置の構成===以下、図1を参照して、本発明の一実施形態における空調システム全体の構成の概略および熱源制御装置の構成について説明する。図1に示されている空調システムは、空調熱源システム2の運転制御を行って、対象となる施設内や室内の温度湿度などの環境を調整するシステムであり、例えば、消費電力の一層の削減が求められているデータセンタなどで用いられる。また、当該空調システムは、空調熱源システム2以外に、熱源制御装置1および温湿度計3を含んで構成されている。空調熱源システム2は、一例として、NR台の冷凍機R1ないしRNR、NC台の冷却塔C1ないしCNC、およびNA台の外気冷房A1ないしANAを備えている。さらに、空調熱源システム2は、搬送機器として、各冷凍機および各冷却塔の共有ポンプPMを備えている。そして、各外気冷房および共有ポンプPMは、冷媒配管などを介して冷却対象負荷Lと接続されている。熱源制御装置1は、計測値データベース11、パラメータテーブル12、時系列予測部13、パラメータ推定部14、混合整数計画演算部15、および制御部16を含んで構成されている。計測値データベース11には、温湿度計3から、時刻k=k0現時刻における外気温度計測値tk0および外気湿度計測値mk0が入力されている。また、計測値データベース11には、現時刻k0における空調熱源システム2の空調負荷計測値yk0も入力されている。なお、計測値データベース11には、所定の計測周期で計測されたこれらの計測値が時系列データとして格納され、当該計測値の時系列データは、時系列予測部13に適宜読み出されている。時系列予測部13からは、所定期間における外気温度、外気湿度、および空調負荷の予測値が時系列データとして出力されている。ここで、所定期間は、計測周期を1ステップとして、現時刻の次のステップk=k01から所定のステップ数先のステップまでの期間であり、当該所定のステップ数をKとすると、例えば、外気温度の予測値の時系列データは、tk01ないしtk0Kと表される。そして、これらの予測値の時系列データは、パラメータ推定部14に入力されている。パラメータテーブル12には、外気温度および外気湿度と対応付けられた空調熱源システム2のモデルパラメータが予め格納されており、当該パラメータテーブル12は、パラメータ推定部14に適宜参照されている。また、後述するように、パラメータ推定部14からは、パラメータ推定部14によって推定された空調熱源システム2のモデルパラメータが反映されたベクトルfk0、hk0、および行列Gk0が出力されている。そして、これらの行列およびベクトルは、混合整数計画演算部15に入力されている。混合整数計画演算部15からは、混合整数計画演算部15によって求められた、所定の制御周期ごとの熱源機器の運転状態を表す状態変数が出力されている。後述するように、本実施形態では、制御周期を計測周期と同一周期とし、混合整数計画演算部15からは、現時刻の次のステップにおける状態変数ベクトルuk01およびδk01が出力される。また、これらの状態変数は、制御部16に入力され、制御部16からは、空調熱源システム2の熱源機器冷凍機、冷却塔、外気冷房および搬送機器共有ポンプPMの制御信号CNTが出力されている。===混合整数計画問題の定式化===本実施形態の熱源制御装置1は、予め空調熱源システム2の運用コストの最小化問題を混合整数計画問題として定式化しておき、制御周期ごとにこの問題を解くことによって、各熱源機器および搬送機器の運転制御を行うための状態変数を求めている。そこで、まず、一例として図1に示した空調熱源システム2に対して、このような混合整数計画問題を定式化する方法について説明する。目的は空調熱源システム2熱源機器および搬送機器の運用コストの最小化であり、目的関数を時刻k=1,2,...,Kにおける空調熱源システム2の運用コストJとし、以下の式1のように定義する。このとき、時刻kにおける運用コスト特性を、以下の式2のように、時変パラメータを有する線形離散事象モデルで表す。ここで、冷凍機、冷却塔、および外気冷房の出力冷却能力を、それぞれの連続変数で表し、時刻kにおける制御入力ベクトルukを、以下の式4のように定義する。また、冷凍機、冷却塔、外気冷房、および共有ポンプPMの起動/停止1/0を、それぞれの0-1変数で表し、時刻kにおける起動/停止ベクトルδkを、以下の式6のように定義する。なお、図1に示した空調熱源システム2においては、共有ポンプの台数NP=1となっている。次に、制約条件として、時刻kにおける空調負荷yk∈Rnと、冷却能力ukおよび起動/停止δkとは、各時刻kにおいて需給バランスを満たす必要があるため、以下の制約条件式7を満たす。また、各熱源機器の冷却能力は以下の上下限制約式8を満たす。さらに、熱源機器の急激な起動/停止の切り替えを防ぐため、以下の連続起動時間制約式9および以下の連続停止時間制約式10を課す。一方、搬送機器である共有ポンプPMは、それを利用する熱源機器の何れかが起動する際、同時に起動される必要があるため、以下の運用制約式11を満たす。ここで、共有ポンプPMを利用する熱源機器の冷却能力の和の最大値をと定義している。ただし、太字のekは太字のukと同じサイズの単位ベクトルである。以上を整理すると、空調熱源システム2の最適制御問題、すなわち、運用コストの最小化問題は、以下の式13ないし式15のように、混合0-1整数計画問題として定式化される。ここで、とする。また、とすれば、である。なお、式18において、diagはdiagonalmatrixの略号であり、対角成分のみ保持する対角行列を意味する。また、式16と式18における、斜体の0は零行列であり、該当する箇所行と列成分が全て0により構成される。===熱源制御装置の動作の概略===以下、図2を適宜参照して、本実施形態における熱源制御装置の動作の概略について説明する。計測値データベース11には、温湿度計3から現時刻k0における外気温度tk0および外気湿度mk0が入力され、空調熱源システム2から現時刻k0における空調負荷yk0が入力される。そして、計測値データベース11は、これらの計測値を所定の計測周期サンプリング周期の時系列データとして格納する。前述したように、パラメータテーブル12には、空調熱源システム2のモデルパラメータが予め格納されている。より具体的には、式13ないし式15のように定式化された混合整数計画問題では、パラメータテーブル12は、時変パラメータa,b,c,dを外気温度tおよび外気湿度mと対応付けた外気温度湿度テーブルとなっている。熱源制御装置1のうち、時系列予測部13、パラメータ推定部14、および混合整数計画演算部15の機能は、例えば、計測値データベース11およびパラメータテーブル12を格納する記憶部と、制御部16を備えるコンピュータによって実現することができる。図2は、時系列予測部13、パラメータ推定部14、および混合整数計画演算部15に相当する機能をコンピュータに実現させるためのプログラムの動作の概略を示している。プログラムの処理が開始されるとS1、まず、記憶部に格納された計測値データベース11を参照して、前述した所定期間時刻k01から時刻k0Kまでにおける外気温度、外気湿度、および空調負荷を予測し、当該Kステップ分の予測値を時系列データとして求めるS2。したがって、当該S2の処理は、時系列予測部13に相当する機能時系列予測処理を実現する。時系列予測処理においては、例えば、公知の記憶ベース推論を用いることができる。当該記憶ベース推論に基づく予測方法では、例えば、現時刻k0からK’ステップ前までの計測値列直近計測値列を計測値データベース11から抽出し、さらに、直近計測値列と類似する同じ時間帯の計測値列類似計測値列を、最小二乗法などを用いて抽出する。そして、抽出した類似計測値列の次のKステップ分の計測値列を直近計測値列と連続するように結合することによって、予測計測値列予測値の時系列データを求めることができる。なお、K’=Kとしてもよい。次に、記憶部に格納されたパラメータテーブル12を参照して、空調熱源システム2のモデルパラメータを推定するS3。したがって、当該S3の処理は、パラメータ推定部14に相当する機能パラメータ推定処理を実現する。より具体的には、時系列予測処理によって予測された時刻k=k01,...,k0Kにおける外気温度予測値tkおよび外気湿度予測値mkに対応する外気温度湿度テーブルを選択し、混合整数計画問題の時変パラメータa,b,c,dを決定する。そして、これらの時変パラメータを式16および式18に当てはめることによって、推定されたモデルパラメータが反映された、時変パラメータベクトルfk0および行列Gk0を求めることができる。一方、時変パラメータベクトルhk0は、時刻kにおける空調負荷予測値ykを当てはめることによって求めることができる。最後に、パラメータ推定処理によって時変パラメータが決定された混合整数計画問題を解き、少なくとも、現時刻の次のステップにおける状態変数ベクトルuk01およびδk01を決定しS4、処理を終了するS5。したがって、当該S4の処理は、混合整数計画演算部15に相当する機能混合整数計画演算処理を実現する。なお、混合整数計画演算処理についての詳細な説明は後述する。制御部16は、混合整数計画演算処理によって求められた状態変数に基づいて制御信号CNTを出力し、空調熱源システム2の各熱源機器および搬送機器の運転制御を行う。そして、以上の処理を制御周期ごとに繰り返すことによって、熱源制御装置1は、空調熱源システム2の運用コストを最小化するように各機器の運転制御を行うことができる。===混合整数計画演算部の動作===以下、図3を参照して、混合整数計画演算部15の動作混合整数計画演算処理について説明する。本実施形態において、混合整数計画演算部15は、図3に示すような分枝限定法を用いて、定式化された混合整数計画問題を解いている。まず、前回の演算結果を利用できない最初の制御周期における動作について説明する。なお、この場合の処理手順は、公知の分枝限定法と同様である。パラメータ推定処理によって混合整数計画問題の時変パラメータが決定されると、混合整数計画演算部15は、当該問題Pをリストに追加するとともに、最適値目的関数値,最適解の暫定値をそれぞれz=∞,x=φ空集合として、混合整数計画演算処理を開始するS401。また、前回の演算結果がなくS402:NO、リストに問題PがあればS404:NO、リストから問題PをPiとして選択するS405。そして、選択された問題Piの緩和問題Pi’を解き、最適解緩和解xi’、およびその場合の最適値z’を求めるS406。緩和問題Pi’を解いた結果、S407ないしS409の判定結果がいずれもNOとなった場合には、状態変数のうち1つの0-1変数を0または1に固定した部分問題Pi1,Pi2を生成し分枝操作、リストに追加する。この際、後述の擬コストを算出するS410。そして、再びS404の判定を行い、リストが空でなければS404:NO、リストの順序に従って新たに問題Piを選択しS405、その緩和問題Pi’を解くS406。一方、緩和問題Pi’を解いた結果、S407ないしS409の何れかの判定結果がYESとなった場合には、問題Piについてこれ以上部分問題を生成して解の探索を行う必要がないため、問題Piをリストから除去し、探索範囲を限定する限定操作。緩和問題Pi’が許容解実行可能解を持たず、実行不可能である場合S407:YESには、問題Piも許容解を持たないため、問題Piをリストから除去しS412、再びS404の判定を行う。また、緩和問題Pi’の最適値z’は、問題Pi、および問題Piについてのすべての部分問題の最適値の下界値となるため、下界値テストと呼ばれるS408の判定を行う。すなわち、z’≧zが成立する場合S408:YESには、問題Piやその部分問題からは、既に得られている最適値の暫定値z≠∞よりも良い小さい最適値は得られないため、問題Piをリストから除去しS412、再びS404の判定を行う。また、緩和問題Pi’の最適解xi’が整数0-1条件を満たし、問題Piの許容解である場合S409:YESには、問題Piの最適解でもあるため、暫定値をそれぞれz=z’,x=xi’として更新するS411。そして、問題Piをリストから除去しS412、再びS404の判定を行う。以上の処理を繰り返し、リストが空になるとS404:YES、その時点の暫定値xを最適解として求め、混合整数計画演算処理を終了するS413。そして、制御部16は、当該最適解xに含まれる状態変数ベクトルuk01およびδk01に基づいて、空調熱源システム2の運転制御を行う。ところで、S410の分枝操作において部分問題を生成する際、0または1に固定する0-1変数分枝変数を何れとするかは、混合整数計画演算処理の効率に大きな影響を与える。本実施形態では、一例として、0-1変数を0または1に固定する優先順位を決定する際に、擬コストと呼ばれる指標値を用いることとする。ここで、j番目の0-1変数の緩和解xj’を0,1に固定した場合の擬コストをそれぞれhdj,hujとすると、目的関数値最適値の増分は、それぞれと推定される。そして、例えば、それぞれの差の絶対値をと求め、djが最大となる0-1変数を優先的に固定して部分問題を生成する。当該部分問題のうち、推定される目的関数値が小さい方の部分問題について先に解の探索を行うことによって、新たな暫定値z=z’が得られた場合S409:YES、S411には、他方の部分問題については、下界値テストによってただちに探索が終了するS408:YES、S412可能性が高くなる。次に、2回目以降の制御周期における動作について説明する。前述したように、計測値データベース11には、所定の計測周期で計測された外気温度湿度などの計測値が入力され、混合整数計画演算部15は、当該計測値を考慮に入れて定式化された混合整数計画問題を、分枝限定法を用いて制御周期ごとに解いている。このように、オンラインで繰り返し分枝限定法を用いて混合整数計画問題を解く場合には、時間的に隣接する制御周期における問題の類似性を活用することができる。本実施形態では、2回目以降の制御周期で、前回直前の制御周期の演算結果がある場合S402:YESには、直前の制御周期において求めた擬コストを現制御周期における擬コストの初期値として設定しS403、解の探索を開始するS404~。したがって、何れの0-1変数を0または1に固定して部分問題を生成すべきかの優先順位情報を初期段階から利用することができるため、解の探索を効率よく行うことができる。また、前回の演算結果を利用する他の方法として、外気温度湿度や空調負荷の変化が微小である場合には、混合整数計画問題の最適解の変化も微小であると想定されるため、新たに混合整数計画問題を解くことなく、前回の演算結果を利用する構成としてもよい。例えば、外気温度湿度や空調負荷の変化量が閾値以下の場合に、制御部16は、直前の制御周期において求めた最適解xに含まれる、時刻k02における状態変数ベクトルuk02およびδk02に基づいて、空調熱源システム2の運転制御を行う構成とすることもできる。このように、前回の演算結果を積極的に利用することによって、混合整数計画演算部15の計算量を低減することができる。前述したように、熱源制御装置1において、予め空調熱源システム2の運用コストの最小化問題を混合整数計画問題として定式化しておき、計測値の時系列データが格納された計測値データベース11を参照して所定期間における外気温度湿度および空調負荷を予測し、これらの予測値に基づいて推定されたモデルパラメータが反映された混合整数計画問題を制御周期ごとに解き、少なくとも、現時刻の次のステップにおける熱源機器の運転状態を表す状態変数を求めることによって、外気冷房を含む熱源機器の外気温度湿度への依存性を考慮した混合整数計画問題を定式化して解くことができる。そのため、混合整数計画法を用いて、外気冷房を備えた空調熱源システム2の運用コストを最小化するように各機器の運転制御を的確に行うことができる。また、混合整数計画問題を解く際に、直前の制御周期における演算結果を利用することによって、計算量を低減することができる。また、各熱源機器の起動/停止を0-1変数で表した混合0-1整数計画問題を分枝限定法を用いて解く際に、擬コストと呼ばれる指標値を求めて、0-1変数を0または1に固定する優先順位を決定するとともに、求めた擬コストを次の制御周期における擬コストの初期値として利用することによって、何れの0-1変数を0または1に固定して部分問題を生成すべきかの優先順位情報を初期段階から利用することができるため、解の探索を効率よく行うことができる。そのため、空調熱源システム2の運用コストの最小化問題について、最適性が保証された厳密な最適解を実用的な計算時間で得ることが可能となる。また、時変パラメータを外気温度湿度と対応付けたパラメータテーブル12を参照することによって、外気温度湿度の予測値に対応する外気温度湿度テーブルを選択し、混合整数計画問題の時変パラメータを決定することができる。また、図1に示した空調システムにおいて、空調熱源システム2および温湿度計3から計測値データベース11に、それぞれ現時刻k0における空調負荷および外気温度湿度の計測値を入力することによって、熱源制御装置1は、当該計測値を考慮に入れて定式化された混合整数計画問題を解き、各熱源機器の運転制御を行うための状態変数を求めることができる。また、熱源機器として、自然エネルギーを活用する外気冷房のほか、冷凍機および冷却塔を組み合わせた空調熱源システム2を構築し、外気温度湿度に応じて各熱源機器の運転制御を行うことによって、外気条件の変動に追従して需給バランスを満たしつつ省エネルギー化を図ることができる。また、熱源制御装置1の時系列予測部13、パラメータ推定部14、および混合整数計画演算部15に相当する機能をコンピュータに実現させるためのプログラムにおいて、予め空調熱源システム2の運用コストの最小化問題を混合整数計画問題として定式化しておき、記憶部に格納された計測値の時系列データに基づいて所定期間における外気温度湿度および空調負荷を予測し、これらの予測値に基づいて推定されたモデルパラメータが反映された混合整数計画問題を制御周期ごとに解き、少なくとも、現時刻の次のステップにおける熱源機器の運転状態を表す状態変数を求めることによって、外気冷房を含む熱源機器の外気温度湿度への依存性を考慮した混合整数計画問題を定式化して解くことができる。そして、制御部16は、求めた状態変数に基づいて、外気冷房を備えた空調熱源システム2の運用コストを最小化するように各機器の運転制御を的確に行うことができる。また、時変パラメータが外気温度湿度と対応付けられたパラメータテーブル12を記憶部に予め格納しておくことによって、外気温度湿度の予測値に対応する外気温度湿度テーブルを選択し、混合整数計画問題の時変パラメータを決定することができる。また、予め空調熱源システム2の運用コストの最小化問題を混合整数計画問題として定式化しておき、計測値の時系列データを参照して所定期間における外気温度湿度および空調負荷を予測し、これらの予測値に基づいて推定されたモデルパラメータが反映された混合整数計画問題を制御周期ごとに解き、少なくとも、現時刻の次のステップにおける熱源機器の運転状態を表す状態変数を求めることによって、外気冷房を含む熱源機器の外気温度湿度への依存性を考慮した混合整数計画問題を定式化して解き、求めた状態変数に基づいて、外気冷房を備えた空調熱源システム2の運用コストを最小化するように各機器の運転制御を的確に行うことができる。なお、上記実施形態は、本発明の理解を容易にするためのものであり、本発明を限定して解釈するためのものではない。本発明は、その趣旨を逸脱することなく、変更、改良され得るとともに、本発明にはその等価物も含まれる。"