(57)【要約】【課題】例えばモデルの変数に対応する誤差項の最大値としての誤差上限値の演算を実用的な時間内や使用メモリ量で実行可能とする。【解決手段】モデルパラメータ決定プログラムは、変数およびモデルパラメータの初期値を用いて計算された変数、変数の微分値をモデルに対応する一階述語論理式に代入する手順と、代入後の論理式内の未知パラメータに仮定値を代入する手順と、その後の論理式に限定記号消去法を適用して誤差上限値を求める手順とを計算機に実行させる。【選択図】図１
2007272478.tif 000002
【特許請求の範囲】
【請求項１】
モデルパラメータの値を決定する計算機によって使用されるプログラムであって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入する手順と、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入する手順と、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求める手順とを計算機に実行させることを特徴とするモデルパラメータ決定プログラム。
【請求項２】
前記モデルパラメータの仮定値代入手順において、前記モデルの変数のうちで、少なくとも１つ以上を除く変数に対する近似値を、さらに前記一階述語論理式に代入することを特徴とする請求項１記載のモデルパラメータ決定プログラム。
【請求項３】
前記求められた誤差上限値を初期値として、該誤差上限値を最小値に絞り込む手順をさらに計算機に実行させることを特徴とする請求項１記載のモデルパラメータ決定プログラム。
【請求項４】
前記最小値への絞込み手順において求められた誤差上限値の最小値を用いて、限定記号消去法を適用して該一階述語論理式に対応するモデル内の前記未知モデルパラメータの値を決定する手順をさらに計算機に実行させることを特徴とする請求項３記載のモデルパラメータ決定プログラム。
【請求項５】
モデルパラメータの値を決定する装置であって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入する手段と、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入する手段と、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求める手段とを備えることを特徴とするモデルパラメータ決定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、理工学の各種の分野において一階述語論理式、すなわち限定記号付の式によって表現される広範な範囲の応用問題を取り扱うものであるが、ここでは生体系シミュレーションにおける応用を主要な例として、本発明の内容を説明する。
【０００２】
生体や細胞に関する各種の系、すなわち生体系のモデル、例えば解糖系反応などのモデルに対してシミュレーションを行い、生体系の解析や設計を行う分野において、モデルのパラメータを効率的に決定することが望まれている。
【０００３】
すなわち本発明で扱う１つの問題は、例えばハイブリッドペトリネット（ＨＰＮ）によってモデル化された、対象となる生体系のモデルに対応して数値シミュレーション、または実測により生体系の変数の時系列データの値が与えられたときに、モデルに基づいて生体系の性質、あるいは本質を把握することである。具体的には、例えばモデルのパラメータとしての反応係数の実行可能範囲を求めること、モデルのパラメータのいくつかの実測値がわかっているときに残りのパラメータの実行可能領域を求めること、モデルの解の存在を確認すること、生体系に対する新しい知見を発見することである。
【０００４】
近年バイオ研究が盛んになるにつれて、このような要求を効率的に満たす有効な解決手法の開発が望まれている。そのような手法は、産業上の生体系設計過程、例えば疾病発症メカニズムの解明、予防医療、テーラーメイド医療などで、極めて広範な範囲において必要とされるものである。しかしながら従来においては、このような要求を効率的に満たす有効な手法は存在せず、生体系のモデルを構築して数値シミュレーションを行うためのツール、例えばゲノミックオブジェクトネット（ＧＯＮ）、ビジュアルオブジェクトネット（ＶＯＮ）、Ｅ−ＣＥＬＬなどを用いて、パラメータの値を様々に設定し、数値シミュレーションを繰り返し、トライアルアンドエラーの方式で望ましいパラメータの値を決定したり、系の性質を推測したりする方法が行われており、実際の作業としては大変に困難で、コストのかかる作業となっている。
【０００５】
このようなシミュレーションにおけるモデルパラメータの決定方式に関する従来技術としての特許文献１では、モデルパラメータに関する制約式を生成し、その制約式に対して限定記号消去法、すなわちＱＥ（クウォンティファイヤ・エリミネーション）アルゴリズムを適用して、解の存在の有無や、実行可能なパラメータの存在範囲を求める技術が開示されている。
【０００６】
また非特許文献１では、同様にＱＥアルゴリズムを用いて生化学モデルにおけるパラメータの決定法が開示されている。
【特許文献１】特開２００５−１２９０１５号公報  「シミュレーションにおけるモデルパラメータ決定プログラム、および決定装置」
【非特許文献１】Ｏｒｉｉ，Ｓ．，Ａｎａｉ，Ｈ．ａｎｄ  Ｈｏｒｉｍｏｔｏ，Ｋ．，“Ｓｙｍｂｏｌｉｃ−ｎｕｍｅｒｉｃ  ｅｓｔｉｍａｔｉｏｎ  ｏｆ  ｐａｒａｍｅｔｅｒｓ  ｉｎ  ｂｉｏｃｈｅｍｉｃａｌ  ｍｏｄｅｌｓ  ｂｙ  ｑｕａｎｔｉｆｉｅｒ  ｅｌｉｍｉｎ−ａｔｉｏｎ”，Ｐｒｏｃ．  ｏｆ  ｔｈｅ  ２００５    Ｉｎｔｅｒｎａｔ．  Ｊｏｉｎｔ  Ｃｏｎｆ．  ｏｆ  ＩｎＣｏＢ，ＡＡＳＢｉ  ａｎｄ  ＫＳＢＩ，２７２−２７７，２００５
【０００７】
図１２は、非特許文献１によって開示されているパラメータ決定方式の従来技術のフローチャートである。同図は、次の一階述語論理式によって表わされる系に対するモデルパラメータ決定処理のフローチャートである。
【０００８】
【数１】
2007272478.tif 000003
【０００９】
この一階述語論理式ＦにおいてＸは変数ベクトル、
【００１０】
2007272478.tif 000004
【００１１】
は変数微分値のベクトル、Ｋはパラメータベクトル、ｅｍａｘは本発明で問題とする、例えば各変数に対応する誤差の上限値であり、ｅｉは各変数に対応する誤差項を示す。この誤差項、および誤差の上限値についてはさらに後述する。
【００１２】
またこの一階述語論理式におけるΨは制約条件、ｈｌ（Ｘ）＝０は、例えばエネルギーなどの保存則、Ｘｉ∈Ｄｉ(Ｄｉ=[ａ,ｂ], ａ,ｂ∈Ｒ∪{∞})は変数の可能領域、
【００１３】
【数２】
2007272478.tif 000005
【００１４】
は、測定値
【００１５】
2007272478.tif 000006
【００１６】
とシミュレーションによる計算値ｘｊとが一致すべきことを示す目的関数であり、｜ｅｊ｜≦ｅｍａｘは、例えば各変数に対する誤差項eiの絶対値が誤差の上限値emaxより小さくなければならないことを示す。
【００１７】
図１２において処理が開始されると、まずステップＳ１００で変数ベクトルＸの初期値Ｘ（０）と、パラメータベクトルＫの初期値Ｋ（０）の値が設定され、ステップＳ１０１で数値計算によってベクトルＸと
【００１８】
2007272478.tif 000007
【００１９】
が求められ、ステップＳ１０２で計算された変数ベクトルＸとその微分のベクトル
【００２０】
2007272478.tif 000008
【００２１】
の値が一階述語論理式Ｆに代入され、論理式Ｆ’が得られる。
すなわちステップＳ１００では本来未知であるパラメータベクトルＫの各要素のパラメータに対して初期値が仮定されて、その初期値を用いて数値計算されたベクトルＸと
【００２２】
2007272478.tif 000009
【００２３】
の値が一階述語論理式に代入される。この場合、以後の計算を厳密に行うために、変数ベクトルＸの要素のうちのいくつかの要素については値は代入されず、その他の変数だけに値が代入されて論理式Ｆ’が作られる。すなわちベクトルＸの要素、例えばｘ１からｘｎまでのｎ個のうちで、いくつかの要素から構成されるベクトルＹについては変数の値が代入されないものとする。
【００２４】
例えば測定器のノイズの影響などを考慮するためにオフセット値を導入すべき場合には、ステップＳ１０３でそのオフセットが求められ、ステップＳ１０２で作られた論理式Ｆ’にさらにオフセットが代入されて、Ｆ’として再びＳ１０２の処理に戻される。オフセットを考慮する必要がないときにはステップＳ１０３の処理は実行されない。
【００２５】
続いてステップＳ１０４で論理式Ｆ’にＱＥアルゴリズムが適用されて未知変数ベクトルＹ、パラメータベクトルＫ、および誤差項ｅｉが消去され、誤差の上限値ｅｍａｘの値が求められる。ここでｅｍａｘは、例えばｎ個の変数ｘｉからｘｎまでのそれぞれに対応して、ＱＥアルゴリズムの収束を確実にするために設定される誤差ｅｉ（ｉ＝１，２，．．，ｎ）に対する誤差の上限値（絶対値）を示すものであり、この上限値ｅｍａｘの値はできるだけ小さいことが望ましい。
【００２６】
ｅｍａｘの値が求められると、続いてステップＳ１０５でその値を用いて再びＱＥアルゴリズムが適用され、ステップＳ１００で初期値が仮定された、実際には未知のパラメータが求められ、ステップＳ１０６でそのパラメータの値を用いて各変数の測定値とシミュレーションによる計算値との残差の二乗和が計算され、ステップＳ１０７でその二乗和の値があらかじめ定められた閾値以下であるか否かが判定される。以下でない場合には、ステップＳ１００以降の処理が繰り返され、閾値以下であると判定された場合にはステップＳ１０８でパラメータの値が出力されて処理を終了する。
【００２７】
特許文献１の従来技術でも、例えば各変数に対応して誤差項が導入され、各誤差項の誤差の上限値を示すｅｍａｘが設定されてモデルパラメータの決定が行われる。
一般に数値シミュレーションや実験によるデータを変数の値として用いる場合には、数値計算の誤差や観測誤差が含まれるために、実際には解があっても数式処理によって計算を行うＱＥアルゴリズムにおいては解無しという結果が得られてしまうことがある。そこで従来技術においても、この誤差に対する上限値ｅｍａｘが変数として追加され、変数ｅｍａｘの値がＱＥアルゴリズムによって求められた後に、その値を用いて再びＱＥアルゴリズムが適用されて未知パラメータが決定される。
【００２８】
しかしながらＱＥアルゴリズムを一階述語論理式に適用して誤差の上限値としてのｅｍａｘを求める場合には、変数の数、未知パラメータの数、観測データの数が増えるにしたがって問題の規模が大きくなり、計算時間、および使用メモリ量が増え、ついには解が得られなくなってしまうという問題点がある。すなわちＱＥアルゴリズムにおける演算量は、一般に変数量とともにダブルイクスポーネンシャル的に増加し、問題が大規模になるとともに実用的な計算時間の範囲内で誤差の上限値ｅｍａｘを求めることができなくなるという問題点がある。
【発明の開示】
【発明が解決しようとする課題】
【００２９】
本発明の課題は、上述の問題点に鑑み、ＱＥアルゴリズムを用いて、例えば各変数に対応して設定された誤差項のうちの最大値、すなわち誤差上限値を求める演算を、問題の規模が大きくなっても実用的な時間内に実行することを可能とし、結果的にモデルパラメータ決定処理を効率化することである。
【課題を解決するための手段】
【００３０】
図１は、本発明のモデルパラメータ決定プログラムの原理的な機能ブロック図である。同図において、まずステップＳ１で変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の微分値を、モデルに対応する一階述語論理式に代入する手順が実行され、ステップＳ２でその一階述語論理式の内部の未知のモデルパラメータに仮定値を代入する手順が実行され、ステップＳ３でその一階述語論理式に限定記号消去法が適用され、その一階述語論理式内の複数の誤差項に共通する誤差上限値ｅｍａｘを求める手順が計算機によって実行される。
【００３１】
発明の実施の形態においては、ステップＳ２の未知パラメータへの仮定値代入手順において、ステップＳ１で、例えばｎ個の変数のうちで、前述のベクトルＹの要素、すなわち値が代入されなかった変数のうちで、少なくとも１つ以上を除く変数に対してさらに値を代入し、実質的に変数の値をさらに減少させて、誤差上限値ｅｍａｘの計算をさらに効率化することも可能である。
【００３２】
また実施の形態においては、ステップＳ３で求められた誤差上限値を初期値として誤差上限値を最小値に絞り込む手順と、その絞込み手順において求められた誤差上限値の最小値を用いて、限定記号消去法により未知パラメータの値を決定する手順をさらに計算機に実行させることも可能である。
【００３３】
次に本発明のモデルパラメータ決定装置は、図１のステップＳ１の処理に相当する動作を行う手段と、ステップＳ２の処理に相当する動作を行う手段と、ステップＳ３に相当する動作を行う手段とを備える。
【００３４】
以上のように本発明においては、誤差上限値を求めるにあたって、未知パラメータ、および正確な値が不明の変数に対して近似的な仮定値を代入して限定記号消去法が適用される。
【発明の効果】
【００３５】
本発明によれば、問題の規模が大きくなり、変数の数や求めるパラメータの数が増大しても、誤差上限値の値を実用的な時間内に近似的に求めることが可能となり、さらにその誤差上限値を初期値として最小値への絞込みを行って、絞り込まれた最小値を用いて未知パラメータの値を決定することが可能となり、モデルパラメータ決定処理の効率化に寄与するところが大きい。
【発明を実施するための最良の形態】
【００３６】
本発明においては、例えばバイオ生化学反応のメカニズムの解明のためにモデルパラメータの決定を行うに当たって、例えば一階述語論理式内で各変数に対して設定される誤差項の上限値を実用的な時間内で近似的に限定記号消去法を用いて求めることによって、モデルパラメータの決定処理を効率化することが可能となるものであり、まずこの限定記号消去法についてその概要を説明する。
【００３７】
多くの産業上の問題や数学の問題は方程式、不等式、限定記号、ブール演算記号などからなる式として記述される。このような式は一階述語論理式と呼ばれ、限定記号消去法（ＱＥ）のアルゴリズムは、与えられた一階述語論理式に対して、限定記号のない等価な式を構成するためのアルゴリズムである。
【００３８】
この限定記号消去（Ｑｕａｎｔｉｆｉｅｒ  Ｅｌｉｍｉｎａｔｉｏｎ）法についてはその概要を紹介した次の文献がある。
【非特許文献２】穴井宏和「Ｑｕａｎｔｉｆｉｅｒ  Ｅｌｉｍｉｎａｔｉｏｎ−アルゴリズム・実装・応用−」  数式処理  Ｖｏ１．１０，Ｎｏ．１，ｐｐ．３−１２（２００３）
【００３９】
図２はこのＱＥ法の概要の説明図である。同図において入力は多項式や不等式などを用いた一階述語論理式であり、出力は限定記号のないパラメータの実行可能領域であり、全ての変数が限定されている場合にはそのｔｒｕｅ（真）、またはｆａｌｓｅ（偽）、すなわち解が存在するかしないか、また存在する場合にはサンプルの解も出力として得ることができる。このような問題は決定問題といわれる。
【００４０】
図２において変数ｘについてｘ２＋ｂｘ＋ｃ＝０という限定のついた問題に対して、限定記号のない等価な式としてｂ２−４ｃ≧０という式が得られている。
いくつかの変数に対して限定記号がない場合には、ＱＥアルゴリズムによって与えられた一階述語論理式と等価な限定記号のない式が得られるが、得られた式は限定記号がついていない残りの変数の可能な範囲を記述するものである。そのような範囲が存在しない場合、ｆａｌｓｅが出力される。このような問題は一般限定記号消去法問題と呼ばれる。
【００４１】
図３は、この限定記号消去法の具他的な制約問題の解法への応用例の説明図である。同図において、例１ではｘとｙとの両方について限定記号がついているため、ＱＥアルゴリズムによってｔｒｕｅであることと、サンプル解とが出力されている。例２においては、ｘだけに限定記号がついているため、ＱＥアルゴリズムによってもう１つの変数としてのｙの実行可能範囲が出力されている。
【００４２】
図４は、本実施形態におけるモデルパラメータ決定装置の機能構成ブロック図である。本実施形態においては、例えば生体系のモデル化が生物学の知識を用いて行われ、シミュレーションに適したモデルが得られ、そのモデルに対応してバイオシミュレーションが行われる。
【００４３】
図４においては、計算機システム１０に対して、入力装置１１を介してシミュレーション結果や、後述するパラメータなどの初期値が入力され、また入力装置１３を介して実験値、すなわち観測値が入力され、それぞれファイル１２、１４に格納される。
【００４４】
これらのファイル内のデータは、その値、あるいは実行可能範囲を決定したいパラメータとしての反応係数などに対する定式化を行う定式化部１５に入力される。定式化部１５は、後述するように実験値などに対する誤差を変数として扱うために誤差変数を導入する誤差変数導入部、ある時刻ｔにおける反応の速度を求める時刻ｔ反応速度計算部、制約問題あるいは最適化問題としての制約式を作成する制約式生成部、後述するように変数や誤差変数に対して制約を付け加えるために、決定したい変数と誤差変数を含む式をＥとして抽出するＥ抽出部、およびＥに対する制約付加部を備えている。
【００４５】
ＱＥ部１６、知見取得コンポーネント部１７による処理は、数式処理エンジン部１８を用いて行われる。ＱＥ部１６はｔｒｕｅ（真）、またはｆａｌｓｅ（偽）を出力し、解の有無を判定する解の存在判定部、および反応係数などのパラメータの可能範囲計算部を備え、また知見取得コンポーネント部１７はパラメータの絞込み部と、パラメータ間の関係解析部とを備え、それぞれの処理結果は出力装置１９、例えばプリンタやディスプレイによって出力される。
【００４６】
以下本実施形態における処理について図５、図６のフローチャートを用いて説明する。図５は、図１２で説明した従来技術における誤差上限値の求め方を改良した誤差上限値初期値計算処理のフローチャートである。同図においてステップＳ１０からステップＳ１２（ステップＳ１３の処理を含む）の処理は、従来例の図１２におけると同様に実行される。前述のようにステップＳ１２で一階述語論理式Ｆに代入される変数は、例えばＮ個のｘｉ（ｉ＝１，２，．．，ｎ）のうちで数値計算による値の正確度が高いと考えられる変数に限定され、前述のように正確度が低いと考えられる複数の変数から構成される変数ベクトルＹは未知のものとして残される。
【００４７】
続いてステップＳ１４でモデルパラメータ、本来は未知のモデルパラメータに初期値を代入することによって、ステップＳ１２で求められた一階述語論理式Ｆ’から一階述語論理式Ｆ’’が求められる。
【００４８】
続いてステップＳ１５で、可能であれば未知ベクトルＹの要素としての、ステップＳ１２で値が代入されなかった変数のうちで、近似的な値が得られるものについては、その近似的な値が一階述語論理式Ｆ’’に代入されて、一階述語論理式Ｆ’’’が求められる。このような近似的な値が得られる変数がない場合にはステップＳ１５の処理は実行されない。
【００４９】
最後にステップＳ１６で一階述語論理式Ｆ’’’に対してＱＥアルゴリズムが適用され、誤差上限値ｅｍａｘの初期値が求められて処理を終了する。従来技術を示す図１２では、図５におけるステップＳ１４、Ｓ１５の処理を行うことなく、ステップＳ１０４で一階述語論理式Ｆ’に対してＱＥアルゴリズムが適用され、誤差上限値ｅｍａｘを求める処理が実行されるが、前述のように問題の規模が大きくなり、パラメータの数や変数の数が増大すると、計算時間、および使用メモリ量が爆発的に増加し、ついには解を得られなくなるという問題点があったが、本実施形態においては、ステップＳ１４で未知パラメータに対して仮定された値が代入され、さらに可能であればステップＳ１５で未知ベクトルＹの要素の変数のうちで、近似的な値を代入可能な変数については、その値が代入された後にＱＥアルゴリズムが適用されるために、近似的ではあっても実用的な計算時間の範囲内で誤差上限値ｅｍａｘを求めることが可能となる。そしてこのｅｍａｘの値を初期値ｅｍａｘ（０）として、次の図６において誤差上限値の最小値への絞込みを行うことによって、誤差上限値の最小値を求めることができる。
【００５０】
図６は、その誤差上限値の最小値計算処理の詳細フローチャートである。同図において処理が開始されると、まずステップＳ２０で図５のステップＳ１６で求められた誤差上限値ｅｍａｘ（０）を、例えば整数ａで割った値をｅｍａｘ、また新たな変数（ｅｍａｘ）Ｌを０とする初期化処理が行われた後に、ステップＳ２１以降の処理が行われる。
【００５１】
ステップＳ２１では、ＱＥアルゴリズム適用の対象となる一階述語論理式にｅｍａｘの値が代入され、ステップＳ２２でＱＥアルゴリズムの適用が行われ、ステップＳ２３でその適用結果が真であるか、偽であるかが判定される。
【００５２】
ＱＥアルゴリズムの適用結果が真でない場合には、ステップＳ２7でｅｍａｘを（ｅｍａｘ）Ｌとし、ｅｍａｘとｅｍａｘ（０）の中間の値をｅｍａｘとする処理が行われた後に、ステップＳ２１以降の処理が続行される。
【００５３】
ステップＳ２３でＱＥアルゴリズムの適用結果が真である場合には、ステップＳ２４で得られたｅｍａｘに対する相対誤差があらかじめ定められた閾値θ未満になっているか否かが判定される。ここでは、この判定は、例えば次式に従って行われる。
【００５４】
【数３】
2007272478.tif 000010
【００５５】
ステップＳ２４で相対誤差が閾値未満でないと判定されると、ステップＳ２８でｅｍａｘをｅｍａｘ（０）とし、ｅｍａｘと（ｅｍａｘ）Ｌとの中間の値をｅｍａｘとする処理が行われた後に、ステップＳ２１以降の処理が続行され、ステップＳ２４で相対誤差が閾値未満であると判定されると、ｅｍａｘが誤差上限値の最小値として出力され、処理を終了する。なお、図５、図６における計算の具体例については後述する。
【００５６】
以下具体例を用いてモデルパラメータの決定時における誤差上限値の計算についてさらに説明する。図７は、ＨＩＶプロテナーゼの生化学反応のメカニズムの説明図である。同図においてＥはタンパク質分解酵素であり、ＨＩＶを発現させる方向の酵素である。ＩはＨＩＶに対する阻害剤である。ＰはＨＩＶを発現させるウィルスタンパク質であり、Ｓはその前駆体タンパク質（基質）である。Ｅは単独で存在する場合には、Ｍと平衡状態となる。
【００５７】
図７において、例えば１番上の反応式において、右方向の反応速度がｖ１であり、この反応速度は、右方向の速度を決める係数ｋ１１と左方向の反応速度を決める係数ｋ１２とによって決定される。またこれらの酵素などの中でタンパク質分解酵素Ｅと前駆体タンパク質Ｓとは、ｔ＝０で０であると反応が始まらないためにｔ＝０で正の値とされる。阻害剤Ｉは外部から与えられるものであり、ウィルスタンパク質Ｐなどはｔ＝０で０であることもできる。
【００５８】
図７のモデルにおいて、各反応式の速度を決める係数ｋ１１からｋ６までの１０個がモデルパラメータであり、そのフィッティングにおける初期値は、図８で与えられるものとする。だだしこれらのうち、ｋ２２、ｋ３、ｋ４２、ｋ５２、およびｋ６の５つは、ＱＥアルゴリズムにおける決定対象モデルパラメータであり、これらの値を用いることなく、ＱＥアルゴリズムは実行される。
【００５９】
図７における各反応式の反応速度ｖ１からｖ６を用いてモデルは次の常微分方式によって表わされる。
【００６０】
【数４】
2007272478.tif 000011
【００６１】
この微分方程式を差分化し、図８の１０個のパラメータのうちの５個に値を代入した関係式（制約問題における制約式に相当）を以下に示す
Ψ(M,E,S,ES,P,EP,I,EI,EJ,JM,JE,JS,JES,JP,JEP,JI,JEI,JEJ,k22,k3,k42,k52,k6, erm,ere,ers,eres,erp,erep,eri,erei,erej,emax)=
erm + JM + 2*(1/10*M*M-1/10000*E) = 0 and
ere + JE - ((1/10*M*M-1/10000*E)-(100*S*E - k22*ES) + (k3*ES) - (100*E*P - k42*EP) - (100*E*I - k52*EI)) = 0 and
ers + JS + ((100*S*E - k22*ES) = 0 and
eres + JES - ((100*S*E - k22*ES) - (k3*ES)) = 0 and
erp + JP - ((k3*ES) - (100*E*P - k42*EP)) = 0 and 
erep + JEP - (100*E*P - k42*EP) = 0 and
eri + JI + (100*E*I - k52*EI) = 0 and
erei + JEI - ((100*E*I - k52*EI) -k6*EI) = 0 and
erej + JEJ - k6*EI = 0 and
M>=0 and E >0 and S >0 and ES >=0 and P>=0 and EP>=0 and I>0 and EI>=0 and EJ>=0 and  k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
- emax < erm < emax and
- emax < ere < emax and
- emax < ers < emax and
- emax < eres < emax and
- emax < erp < emax and
- emax < erep < emax and
- emax < eri < emax and
- emax < erei < emax and
- emax < erej < emax and
emax >= 0
ここで上から９番目の式までは前述の微分方程式を差分化したものである。一般に数値シミュレーションや実験によるデータを変数の値として用いる場合には、数値計算の誤差や観測誤差が含まれるために、実際には解があっても、数式処理によって正確な計算を行うＱＥアルゴリズムにおいては解無しという結果が得られてしまうことがある。そこで本実施形態では、各差分式の計算値にそれぞれ微小な誤差が含まれているものとし、誤差の変数を用いて差分式を作成し、ＱＥを適用することにする。実際にはそれぞれの差分式に対する誤差変数として、各変数の先頭にerをつけたerm, ere, ers, eres, erp, erep, eri, erei, およびerejの９つを用いて、誤差項を含む差分式を作成した。これらの差分式に、各変数、およびモデルパラメータに対して物理的に自明な範囲を与える不等式と、各誤差項に対してその絶対値の最大値をｅｍａｘとする不等式を追加したものが式として定義される。なお各差分式内のJ*はd*/dtを表わし、例えばJMはdM/dtである。
【００６２】
まず未知パラメータに対して次のように初期値が与えられる。
% Starting parameters
k22 := 300;
k3 := 10;
k42 := 500;
k52 := 1/10;
k6 := 1/10;
続いて図５のステップＳ１３に対応して次のようにオフセットの計算が行われる。
% Compute offset
t3600_msd :=632629/1000000; % =0.632629
ofst := t3600_msd - 3/125*Sinit;
本実施形態においては、観測点（実測値）として、３つの時刻における変数値や変数の微分値などが、ＱＥ入力として図４のＱＥ部１６に与えられるものとする。そしてまず最初の時刻に対応するＱＥ部１６への入力は次のようになる。ここで上から９番目までの式は前述の関係式Ψにおける９番目までの式を表わし、１０番目の式はその後の計算において用いられる変数ｓｇｎｌの値を与えるものである。
【００６３】
erm := - ( JM + 2*(1/10*M*M-1/10000*ee) );
ere := - ( JEE - ((1/10*M*M-1/10000*ee) - (100*S*ee - k22*ES)
+ (k3*ES) - (100*ee*P - k42*EP) - (100*ee*ii - k52*EI)) );
ers := - ( JSS + (100*S*ee - k22*ES) );
eres := - ( JES - ((100*S*ee - k22*ES) - (k3*ES)) );
erp := - ( JPP - ((k3*ES) - (100*ee*P - k42*EP)) );
erep := - ( JEP - (100*ee*P - k42*EP) );
eri := - ( JII + (100*ee*ii - k52*EI) );
erei := - ( JEIEI - ((100*ee*ii - k52*EI) -k6*EI) );
erej := - ( JEJ - k6*EI );
sgnl := 369873/1000000; % data No.=682, t=984, signal=0.369873
uqfin0:=(
-M - 2*ee + 2*S + 2*P + 2*ii = - 2*Einit + 2*Sinit + 2*Iinit and
S + ES + P + EP = Sinit and
M + 2*ee - 2*S - 2*P + 2*EI +2*Ej = 2*Einit - 2*Sinit and
p - (sgnl - ofst)*1000/24 =0 and
M>=0 and ee >0 and S >0 and ES >=0 and P>=0 and EP>=0 and EI>=0 and EJ>=0 and
k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
(- emax <= erm <= emax and
- emax <= ere <= emax and
- emax <= ers <= emax and
- emax <= eres <= emax and
- emax <= erp <= emax and
- emax <= erep <= emax and
- emax <= eri <= emax and
- emax <= erei <= emax and
- emax <= erej <= emax) );
uqfout0 := sub(
time = 9840000000000/1000000000000*10**002, % 9.840000000000e+002
Iinit = 3000000000000/1000000000000/10**003, % 3.000000000000e-003
signal = 1887252952430/1000000000000/10**001, % 1.887252952430e-001
M = 2229085252834/1000000000000/10**005, % 2.229085252834e-005
SS = 1813582723453/1000000000000*10**001, % 1.813582723453e+001
ES = 4876923348145/1000000000000/10**004, % 4.876923348145e-004
PP = 7863553968458/1000000000000*10**000, % 7.863553968458e+000
EP = 1311046795513/1000000000000/10**004, % 1.311046795513e-004
ii = 1330477176374/1000000000000/10**005, % 1.330477176374e-005
EIEI = 5661674744077/1000000000000/10**007, % 5.661674744077e-007
EJ = 2986129060762/1000000000000/10**003, % 2.986129060762e-003
Einit = 3700000000000/1000000000000/10**003, % 3.700000000000e-003
Sinit = 2600000000000/1000000000000*10**001, % 2.600000000000e+001
Jtime = 5000000000000/(10**13)*(2**001), % 1.0000000000000
JIinit = 0, % 0.0000000000000
Jsignal = 9571837680298/(10**13)/(2**013), % 0.0001168437217
JM = 5560137236768/(10**13)/(2**025), % 0.0000000165705
Jee = -7348462184181/(10**13)/(2**029), % -0.0000000013688
JSS = -6231588486400/(10**13)/(2**007), % -0.0048684285050
JES = -5825875387525/(10**13)/(2**022), % -0.0000001388997
JPP = 6231665156586/(10**13)/(2**007), % 0.0048684884036
JEP = 6627088471687/(10**13)/(2**023), % 0.0000000790011
Jii = -8888929747922/(10**13)/(2**024), % -0.0000000529822
JEIEI = -6078516796086/(10**13)/(2**028), % -0.0000000022644
JEJ = 9268837038598/(10**13)/(2**024), % 0.0000000552466
JEinit = 0, % 0.0000000000000
JSinit = 0, % 0.0000000000000
uqfin0 );
【００６４】
続いて次の時刻に対応するQE部１６への入力は次のようになる。
erm := - ( JM + 2*(1/10*M*M-1/10000*ee1) );
ere := - ( JEE - ((1/10*M*M-1/10000*ee1) - (100*s1*ee1 - k22*ES)
+ (k3*ES) - (100*ee1*p1 - k42*EP) - (100*ee1*ii - k52*ei1)) );
ers := - ( JSS + (100*s1*ee1 - k22*ES) );
eres := - ( JES - ((100*s1*ee1 - k22*ES) - (k3*ES)) );
erp := - ( JPP - ((k3*ES) - (100*ee1*p1 - k42*EP)) );
erep := - ( JEP - (100*ee1*p1 - k42*EP) );
eri := - ( JII + (100*ee1*ii - k52*ei1) );
erei := - ( JEIEI - ((100*ee1*ii - k52*ei1) -k6*ei1) );
erej := - ( JEJ - k6*ei1 );
sgnl := 36621/250000; % data (1) No.=628, t=336, signal=0.146484
uqfin1:=(
-M - 2*ee1 + 2*s1 + 2*p1 + 2*ii = - 2*Einit + 2*Sinit + 2*Iinit and
s1 + ES + p1 + EP = Sinit and
M + 2*ee1 - 2*s1 - 2*p1 + 2*ei1 +2*Ej = 2*Einit - 2*Sinit and
p1 - (sgnl - ofst)*1000/24 =0 and
M>=0 and ee1 >0 and s1 >0 and ES >=0 and p1>=0 and EP>=0 and ei1>=0 and EJ>=0 and
k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
(- emax <= erm <= emax and
- emax <= ere <= emax and
- emax <= ers <= emax and
- emax <= eres <= emax and
- emax <= erp <= emax and
- emax <= erep <= emax and
- emax <= eri <= emax and
- emax <= erei <= emax and
- emax <= erej <= emax) );
uqfout1 := sub(
time = 3360000000000/1000000000000*10**002, % 3.360000000000e+002
Iinit = 3000000000000/1000000000000/10**003, % 3.000000000000e-003
signal = 1007025285443/1000000000000/10**001, % 1.007025285443e-001
M = 1090949476872/1000000000000/10**005, % 1.090949476872e-005
SS = 2180325350548/1000000000000*10**001, % 2.180325350548e+001
ES = 7216956684414/1000000000000/10**004, % 7.216956684414e-004
PP = 4195938689347/1000000000000*10**000, % 4.195938689347e+000
EP = 8610950859596/1000000000000/10**005, % 8.610950859596e-005
ii = 2158705895960/1000000000000/10**004, % 2.158705895960e-004
EIEI = 1142462779546/1000000000000/10**005, % 1.142462779546e-005
EJ = 2772704782609/1000000000000/10**003, % 2.772704782609e-003
Einit = 3700000000000/1000000000000/10**003, % 3.700000000000e-003
Sinit = 2600000000000/1000000000000*10**001, % 2.600000000000e+001
Jtime = 5000000000000/(10**13)*(2**001), % 1.0000000000000
JIinit = 0, % 0.0000000000000
Jsignal = 7035347335850/(10**13)/(2**012), % 0.0001717614096
JM = 6833982121443/(10**13)/(2**025), % 0.0000000203669
Jee = -9043060580504/(10**13)/(2**023), % -0.0000001078017
JSS = -9159408459733/(10**13)/(2**007), % -0.0071557878592
JES = -5201440220687/(10**13)/(2**019), % -0.0000009920960
JPP = 9160608510293/(10**13)/(2**007), % 0.0071567253987
JEP = 9153008209721/(10**13)/(2**024), % 0.0000000545562
Jii = -5427209515084/(10**13)/(2**019), % -0.0000010351581
JEIEI = -5593339296621/(10**13)/(2**023), % -0.0000000666778
JEJ = 5776793220764/(10**13)/(2**019), % 0.0000011018359
JEinit = 0, % 0.0000000000000
JSinit = 0, % 0.0000000000000
uqfin1 );
【００６５】
さらに最後の３つ目の時刻に対応するＱＥ部１６への入力は次のようになる。
erm := - ( JM + 2*(1/10*M*M-1/10000*ee2) );
ere := - ( JEE - ((1/10*M*M-1/10000*ee2) - (100*s2*ee2 - k22*ES)
+ (k3*ES) - (100*ee2*p2 - k42*EP) - (100*ee2*ii - k52*ei2)) );
ers := - ( JSS + (100*s2*ee2 - k22*ES) );
eres := - ( JES - ((100*s2*ee2 - k22*ES) - (k3*ES)) );
erp := - ( JPP - ((k3*ES) - (100*ee2*p2 - k42*EP)) );
erep := - ( JEP - (100*ee2*p2 - k42*EP) );
eri := - ( JII + (100*ee2*ii - k52*ei2) );
erei := - ( JEIEI - ((100*ee2*ii - k52*ei2) -k6*ei2) );
erej := - ( JEJ - k6*ei2 );
sgnl := 566101/1000000; % data No.=754, t=1848, signal=0.566101
uqfin2:=(
-M - 2*ee2 + 2*s2 + 2*p2 + 2*ii = - 2*Einit + 2*Sinit + 2*Iinit and
s2 + ES + p2 + EP = Sinit and
M + 2*ee2 - 2*s2 - 2*p2 + 2*ei2 +2*Ej = 2*Einit - 2*Sinit and
p2 - (sgnl - ofst)*1000/24 =0 and
M>=0 and ee2 >0 and s2 >0 and ES >=0 and p2>=0 and EP>=0 and ei2>=0 and EJ>=0 and
k22>0 and k3>0 and k42>0 and k52>0 and k6>0 and
(- emax <= erm <= emax and
- emax <= ere <= emax and
- emax <= ers <= emax and
- emax <= eres <= emax and
- emax <= erp <= emax and
- emax <= erep <= emax and
- emax <= eri <= emax and
- emax <= erei <= emax and
- emax <= erej <= emax) );
uqfout2 := sub(
time = 1848000000000/1000000000000*10**003, % 1.848000000000e+003
Iinit = 3000000000000/1000000000000/10**003, % 3.000000000000e-003
signal = 2797735529326/1000000000000/10**001, % 2.797735529326e-001
M = 3668866233730/1000000000000/10**005, % 3.668866233730e-005
SS = 1434217229903/1000000000000*10**001, % 1.434217229903e+001
ES = 3965132539191/1000000000000/10**004, % 3.965132539191e-004
PP = 1165723137219/1000000000000*10**001, % 1.165723137219e+001
EP = 1998155229046/1000000000000/10**004, % 1.998155229046e-004
ii = 3777106832261/1000000000000/10**007, % 3.777106832261e-007
EIEI = 1652850731302/1000000000000/10**008, % 1.652850731302e-008
EJ = 2999605760809/1000000000000/10**003, % 2.999605760809e-003
Einit = 3700000000000/1000000000000/10**003, % 3.700000000000e-003
Sinit = 2600000000000/1000000000000*10**001, % 2.600000000000e+001
Jtime = 5000000000000/(10**13)*(2**001), % 1.0000000000000
JIinit = 0, % 0.0000000000000
Jsignal = 7784861579264/(10**13)/(2**013), % 0.0000950300486
JM = 5662307413375/(10**13)/(2**025), % 0.0000000168750
Jee = 5336027763135/(10**13)/(2**027), % 0.0000000039757
JSS = -5068251389866/(10**13)/(2**007), % -0.0039595713983
JES = -7645065759424/(10**13)/(2**023), % -0.0000000911363
JPP = 5068269257599/(10**13)/(2**007), % 0.0039595853575
JEP = 6474129079970/(10**13)/(2**023), % 0.0000000771776
Jii = -8297377108819/(10**13)/(2**029), % -0.0000000015455
JEIEI = -5745544012119/(10**13)/(2**033), % -0.0000000000669
JEJ = 8656473925608/(10**13)/(2**029), % 0.0000000016124
JEinit = 0, % 0.0000000000000
JSinit = 0, % 0.0000000000000
uqfin2 );
【００６６】
続いてこのようなＱＥ部１６への入力に対するＱＥアルゴリズムの適用結果が次のようにチェックされる。
%Check .true. or .false.
a:=ex({ee,ee1,ee2,S,P,EI,S1,P1,EI1,S2,P2,EI2,emax},uqfout0 and uqfout1 and uqfout2)$
f_a:=rlqe(a);
%Find emax
a1:=ex({ee,ee1,ee2,S,P,EI,S1,P1,EI1,S2,P2,EI2},uqfout0 and uqfout1 and uqfout2)$
f_a1:=rlqe(a1);
emax:=(-part(f_a1,1,2)/part(f_a1,1,1,1));
%Check .true. or .false.
a2:=ex({ee,ee1,ee2,S,P,EI,S1,P1,EI1,S2,P2,EI2},uqfout0 and uqfout1 and uqfout2)$
f_a2:=rlqe(a2);
end$
【００６７】
すなわち、まず第１にＱＥアルゴリズムの適用結果が真であるか偽であるかが判定され、次に誤差上限値ｅｍａｘが求められる。ここでは最後に従来例の図１２におけるステップＳ１０５に対応して、計算された誤差上限値ｅｍａｘの値を用いて再びＱＥアルゴリズムが適用され、未知パラメータを求める処理が行われ、その結果が真になるが偽になるかがチェックされているが、本実施形態においては、その直前で計算された誤差上限値ｅｍａｘの値が、図５のステップＳ１６で求められたｅｍａｘの初期値ｅｍａｘ（０）に当たり、この初期値を用いて図６の処理フローチャートに従って誤差上限値ｅｍａｘの最小値への絞込みが行われる。
【００６８】
次に図５、図６の処理フローチャートに関連して、誤差上限値ｅｍａｘの初期値の計算時間短縮効果と、初期値が求められた後の誤差上限値ｅｍａｘの最小値への絞込みの計算時間について、図９と図１０を用いて説明する。図９は、誤差上限値の初期値の計算時間に対する短縮効果の説明図である。同図においてＮｏ．１は、従来例の図１２におけると同様の処理における計算時間である。ここでは、観測点の数ｎｄ＝３、すなわち３つの時刻における観測値、値が代入されなかった変数の数、すなわちベクトルＹの要素の数ｎＹ＝３、未知パラメータの数ｎＫ＝５に対応する計算時間であり、その値は８０．５秒となっている。
【００６９】
これに対してＮｏ．２は、５個の未知パラメータの全てに対して仮定された値を代入した場合の計算時間であり、未知パラメータの数ｎＫを０とすることによって計算時間はＮｏ．１の１／２４５となっている。
【００７０】
さらにＮｏ．３は、５個の未知パラメータへの値の代入に加えて、未知ベクトルＹの要素としての３つの変数のうちで、２つの変数に近似値を代入した場合の時間であり、Ｎｏ．１の１／２５０００の時間となっている。
【００７１】
このように観測点と数ｎｄと、値が代入されなかった変数の数ｎＹと、未知パラメータの数ｎＫによって次式で求められるｎの値がＮｏ．１では１４、Ｎｏ．２では９、Ｎｏ．３では３と小さくなるに従って、十分実用的な計算時間の範囲で誤差上限値ｅｍａｘの初期値が得られることがわかる。
【００７２】
ｎ＝ｎＫ＋ｎＹ×ｎｄ
ただし、ｅｍａｘの初期値を求めるためのＱＥアルゴリズムの演算における消去対象として、例えば未知変数ベクトルＹの要素のうちで少なくとも1つ以上の変数を、値が代入されないままに残しておく必要がある。
【００７３】
一般的には図９で説明したようにパラメータに対して仮定される値、未知変数に対して近似値を代入することによって得られた誤差上限値ｅｍａｘの値は、従来技術によって求められた値より大きくなる。まれに等しくなることもあるが、一般的には従来技術で求められた値、すなわち最小値に近づける絞込み処理を実行する必要がある。この絞込みの方法としては、各種の方法が考えられるが、本実施形態では二分探索法を用いることとする。
【００７４】
図１０は、図６で説明した二分探索処理による計算時間の説明図である。まずＮｏ．１は、図９において説明した本実施形態において得られた誤差上限値ｅｍａｘの初期値ｅｍａｘ（０）を０．０９９２とし、図６のステップＳ２０におけるａの値を１０として、ステップＳ２２においてＱＥアルゴリズムを適用した結果であり、計算時間は０．１８秒であるが、その結果は偽である。
【００７５】
Ｎｏ．２では、図６のステップＳ２７で０．０９９２と０．００９９２の中間の値の０．０５４５６をｅｍａｘ（０）とし、ＱＥアルゴリズムを適用した結果であり、計算時間は０．２秒であるが、適用結果はやはり偽である。
【００７６】
Ｎｏ．３では、さらにステップＳ２７で０．０９９２と０．０５４５６の中間の値０．０７６８８をｅｍａｘの値としてＱＥアルゴリズムを適用した結果であり、適用結果は真となっている。図６のステップＳ２３で適用結果が真と判定されても、ステップＳ２４で相対誤差が閾値未満でないと判定されると、ステップＳ２８の処理が行われた後に、ステップＳ２１以降の処理が繰り返される。ここではｅｍａｘの値が小さい方に振られ、ＱＥアルゴリズムが適用されるために、ステップＳ２３でその適用結果が再び偽であると判定されることもあるが、いずれにしてもステップＳ２１からＳ２８の処理が繰り返されるうちに、ステップＳ２４で相対誤差が閾値未満と判定された時点で、ステップＳ２５で誤差上限値の最小値が出力され、ステップＳ２６でその最小値を用いてモデルパラメータが決定される。
【００７７】
図１０で説明したＮｏ．３までの計算においても、Ｎｏ．２とＮｏ．３の計算結果で小数点以下第２位の同じ位置に最初の有効数字が現れており、有効数字１桁の範囲内にｅｍａｘの値が収束している。そこで図９においてＮｏ．２の計算を用いるとしても、図１０でＮｏ．３までの計算にかかる全体の時間は０．９９９秒となり、有効数字の桁数を小さくすることによって、従来技術による計算時間と比べて誤差上限値の最小値を求めるまでの時間は十分に短くなり、従来技術の方法では誤差上限値の値そのものを得られない場合があることに比較して、実用的な計算時間内で、近似的ではあるが、誤差上限値の最小値を求める有力な方法として本発明の方法を用いることができる。
【００７８】
以上において本発明のモデルパラメータ決定プログラムについてその詳細を説明したが、このプログラムを用いるモデルパラメータ決定装置は、当然一般的なコンピュータシステムを基本として構成することが可能である。図１１はそのようなコンピュータシステム、すなわちハードウェア環境の構成ブロック図である。
【００７９】
図１１においてコンピュータシステムは中央処理装置（ＣＰＵ）２０、リード・オンリ・メモリ（ＲＯＭ）２１、ランダム・アクセス・メモリ（ＲＡＭ）２２、通信インターフェース２３、記憶装置２４、入出力装置２５、可搬型記憶媒体の読取り装置２６、およびこれらの全てが接続されたバス２７によって構成されている。
【００８０】
記憶装置２４としてはハードディスク、磁気ディスクなど様々な形式の記憶装置を使用することができ、このような記憶装置２４、またはＲＯＭ２１に５、図６のフローチャートに示されたプログラムや、本発明の特許請求の範囲の請求項１〜４のプログラムなどが格納され、そのようなプログラムがＣＰＵ２０によって実行されることにより、本実施形態における誤差上限値の初期値決定、最小値への絞込みなどが可能となる。
【００８１】
このようなプログラムは、プログラム提供者２８からネットワーク２９、および通信インターフェース２３を介して、例えば記憶装置２４に格納されることも、また市販され、流通している可搬型記憶媒体３０に格納され、読取り装置２６にセットされて、ＣＰＵ２０によって実行されることも可能である。可搬型記憶媒体３０としてはＣＤ−ＲＯＭ、フレキシブルディスク、光ディスク、光磁気ディスク、ＤＶＤなど様々な形式の記憶媒体を使用することができ、このような記憶媒体に格納されたプログラムが読取り装置２６によって読取られることにより、本実施形態における誤差上限値の最小値を用いるモデルパラメータの決定が可能となる。
【００８２】
（付記１）  モデルパラメータの値を決定する計算機によって使用されるプログラムであって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入する手順と、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入する手順と、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求める手順とを計算機に実行させることを特徴とするモデルパラメータ決定プログラム。（１）
（付記２）  前記モデルパラメータの仮定値代入手順において、前記モデルの変数のうちで、少なくとも１つ以上を除く変数に対する近似値を、さらに前記一階述語論理式に代入することを特徴とする付記１記載のモデルパラメータ決定プログラム。（２）
（付記３）  前記求められた誤差上限値を初期値として、該誤差上限値を最小値に絞り込む手順をさらに計算機に実行させることを特徴とする付記１記載のモデルパラメータ決定プログラム。
（付記４）  前記最小値への絞込み手順において求められた誤差上限値の最小値を用いて、限定記号消去法を適用して該一階述語論理式に対応するモデル内の前記未知モデルパラメータの値を決定する手順をさらに計算機に実行させることを特徴とする付記３記載のモデルパラメータ決定プログラム。（４）
（付記５）  前記最小値への絞込み手順において、前記求められた誤差上限値と、該誤差上限値より小さい値との間を範囲とする二分探索法を用いることを特徴とする付記３記載のモデルパラメータ決定プログラム。
（付記６）  モデルパラメータの値を決定する計算機によって使用される記憶媒体であって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入するステップと、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入するステップと、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求めるステップとを計算機に実行させるモデルパラメータ決定プログラムを格納した計算機読み出し可能可搬型記憶媒体。
（付記７）  モデルパラメータの値を決定する装置であって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入する手段と、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入する手段と、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求める手段とを備えることを特徴とするモデルパラメータ決定装置。（５）
（付記８）  モデルパラメータの値を決定する方法であって、
モデル内の変数、およびモデルパラメータの初期値を用いて計算された変数、および変数の二分値を該モデルに対応する一階述語論理式に代入し、
該代入後の一階述語論理式内の未知モデルパラメータに仮定値を代入し、
該仮定パラメータ値が代入された一階述語論理式に限定記号消去法を適用して、該一階述語論理式内の複数の誤差項に対応する誤差上限値を求めることを特徴とするモデルパラメータ決定方法。
【図面の簡単な説明】
【００８３】
【図１】本発明のモデルパラメータ決定プログラムの原理的な機能ブロック図である。
【図２】限定記号消去（ＱＥ）法の概要の説明図である。
【図３】制約問題に対するＱＥ法の適用例を示す図である。
【図４】モデルパラメータ決定装置の機能構成ブロック図である。
【図５】誤差上限値初期値計算処理の詳細フローチャートである。
【図６】誤差上限値最小値計算処理の詳細フローチャートである。
【図７】ＨＩＶプロテナーゼのメカニズムを説明する図である。
【図８】図７に対するモデルパラメータの値の例を示す図である。
【図９】誤差上限値初期値の計算時間短縮効果を説明する図である。
【図１０】二分探索による誤差上限値最小値計算時間を説明する図である。
【図１１】本発明のプログラムのコンピュータへのローディングを説明する図である。
【図１２】モデルパラメータ決定処理の従来例のフローチャートである。
【符号の説明】
【００８４】
１０  計算機システム
１１  入力装置
１２  シミュレーション結果、パラメータなどの初期値
１３  入力装置
１４  実験値
１５  定式化部
１６  ＱＥ部
１７  知見取得コンポーネント部
１８  数式処理エンジン部
１９  出力装置
２０  中央処理装置（ＣＰＵ）
２１  リード・オンリ・メモリ（ＲＯＭ）
２２  ランダム・アクセス・メモリ（ＲＡＭ）
２３  通信インターフェース
２４  記憶装置
２５  入出力装置
２６  読み取り装置
２７  バス
２８  プログラム提供者
２９  ネットワーク
３０  可搬型記憶媒体
【図１】
2007272478.tif 000012
【図２】
2007272478.tif 000013
【図４】
2007272478.tif 000014
【図５】
2007272478.tif 000015
【図６】
2007272478.tif 000016
【図７】
2007272478.tif 000017
【図８】
2007272478.tif 000018
【図９】
2007272478.tif 000019
【図１０】
2007272478.tif 000020
【図１１】
2007272478.tif 000021
【図１２】
2007272478.tif 000022
【図３】
2007272478.tif 000023