(57)【要約】【課題】一般的に適用できる効率的な最適化の手法を提供する。【解決手段】本発明は、準最適解または最適解を求める問題を数学的にモデル化する高次元領域の関数を最小化するアプローチによって問題を最適化するシステムを提供する。このシステムは、コンピュータ・プログラムを実行するコンピュータ・システムと、最適化プログラムと、を有する。この最適化プログラムは、最適化問題を受け取り、連続型反復パラメータにより前記最適化問題をパラメータ化して連続最適化問題を生成し、連続最適化問題の解法を採用して前記最適化問題の最適解に向かう軌跡を定め、最適化問題の複数のパラメータを制御変数として定式化し、制御変数を計算するための反復法により最適な軌跡を含む所定領域内で準最適軌跡を求め、準最適軌跡を繰り返し計算して大域的な最適解を含む所定区間内の準最適解を得る。【選択図】図２【特許請求の範囲】
【請求項１】
最適化システムであって、
コンピュータ・プログラムを実行するコンピュータ・システムと、
最適化プログラムと、を有し、
前記最適化プログラムは、
最適化問題を受け取り、
連続型反復パラメータにより前記最適化問題をパラメータ化して連続最適化問題を生成し、
前記連続最適化問題の解法を採用して前記最適化問題の最適解に向かう軌跡を定め、
前記最適化問題の複数のパラメータを制御変数として定式化し、前記制御変数を計算するための反復法により最適な軌跡を含む所定領域内で準最適軌跡を求め、
前記準最適軌跡を繰り返し計算して大域的な最適解を含む所定区間内の準最適解を得る、
最適化システム。
【請求項２】
前記最適化プログラムが、記号形式で表現された最適化問題を受け取る、請求項１に記載の最適化システム。
【請求項３】
前記最適化プログラムが、前記記号形式で表現された最適化問題を受け取った後で、
前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
【請求項４】
前記最適化プログラムが、
前記受け取った最適化問題内の非実数の決定変数を実数の決定変数に変換することによって、前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
【請求項５】
前記受け取った最適化問題内の非実数の決定変数を実数の決定変数に変換することが、
浮動小数点変数を強制的に２つの値０または１のいずれかにするための追加の制約ｙｉ（１−ｙｉ）＝０を伴って範囲［０，１］において、ブール型の各決定変数ｂｉを等価な浮動小数点値ｙｉに変換することと、
整数型の各決定変数を、等価なｋ項の２進多項式に変換することを含み、ここでｋは前記整数型の決定変数がとり得る最大値の対数（基底は２）にほぼ同等であり、
さらに、前記２進多項式をｋ個の浮動小数点変数ｙ０，ｙ１，．．．，ｙｋ−１によるｋ個の項で表すことを含み、前記変数のすべてが範囲［０，１］に含まれ、ｙｉ（１−ｙｉ）＝０の形のｋ個の制約を伴う、
請求項４に記載の最適化システム。
【請求項６】
前記最適化プログラムが、
実数型の決定変数を浮動小数点型の決定変数に変換することによって、前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
【請求項７】
前記最適化プログラムが、
等式制約を不等式制約に変換することによって前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
【請求項８】
前記等式制約を不等式制約に変換することが、
個々の等式制約ｈｊ（ｘ）について、
【数１】
2007509437.tif 000172
のように、ｈｊ（ｘ）を２つの不等式制約に変換し、第３の制約ｓｊ≧０を追加することをさらに含み、ここで、ｍは等式制約を不等式制約に変換する前の不等式制約の数、ｓｊは束縛変数である、請求項７に記載の最適化システム。
【請求項９】
前記最適化プログラムが、
エンベロープ変数を追加すること、および
前記最適化問題の最適解に向かう軌跡への制約の影響を表すバリア関数を定式化することによって、
前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
【請求項１０】
前記エンベロープ変数を追加することが、
不等式制約を伴うエンベロープ変数ｚを、等式制約を不等式制約に変換することによって発生する前記最適化問題
【数２】
2007509437.tif 000173
に加えることをさらに含み、
前記変数ｚは、
【数３】
2007509437.tif 000174
として定義され、前記不等式制約は、
【数４】
2007509437.tif 000175
として定義される、請求項９に記載の最適化システム。
【請求項１１】
前記バリア関数を追加することが、各制約ｇｉについてバリア関数−ｕｉＢ(ｇｉ)を追加することをさらに含む、請求項９に記載の最適化システム。
【請求項１２】
前記バリア関数が、各制約ｇｉについて−ｕｉlog(ｇｉ)の形をとる、請求項１１に記載の最適化システム。
【請求項１３】
前記標準形式による前記最適化問題の記号表現が、
【数５】
2007509437.tif 000176
と表すことができ、ここで、前記関数Ｆは決定変数ｘ、束縛変数ｓ、エンベロープ変数ｚ、複数の制御変数を備えるバリア関数の係数ベクトルｕ、複数の制御変数を備える束縛変数の係数ベクトルｒの関数であり、ｇｉは不等式制約を定義する式であり、
前記最適化問題の前記標準形式による記号表現が、さらにコンパクトに、
【数６】
2007509437.tif 000177
と表すことができ、ここで、
【数７】
2007509437.tif 000178
であり、ここで、ｍは等式制約を不等式制約に変換する前の不等式制約の数、ｋは変換された不等式制約の数、「状態ベクトル」と呼ばれるｙは決定変数ｘ、束縛変数ｓ、およびエンベロープ変数ｚのベクトル、ベクトルｕはバリア関数の係数ベクトル、ｒは束縛変数の係数ベクトルである、請求項３に記載の最適化システム。
【請求項１４】
前記最適化プログラムが、
ｙｋ＋１＝ｙｋ＋φｋ（ｙｋ，ｒｋ，ｕｋ）
（ｋ＝０，１，..．，ｋconvergence）
のように現在の状態ベクトルｙｋから次の状態ベクトルｙｋ＋１を計算できる反復関数φｋをパラメータ化することにより、連続型反復パラメータで前記最適化問題をパラメータ化して連続最適化問題を生成する、請求項１３に記載の最適化システム。
【請求項１５】
前記最適化プログラムが、０から反復限界Ｔまでの範囲にある連続型反復パラメータτを導入し、ベクトル関数ｙ、ｕ、ｒを連続型反復変数τでパラメータ化することによって、前記反復関数φｋをパラメータ化する、請求項１４に記載の最適化システム。
【請求項１６】
前記関数Ｆは、前記連続型反復変数τに関してＦ（ｙ（τ），ｕ（τ），ｒ（τ））と表現される、請求項１５に記載の最適化システム。
【請求項１７】
前記反復関数ｙｋ＋１＝−ｙｋ＋φｋ（ｙｋ，ｒｋ，ｕｋ）が、１次微分方程式
【数８】
2007509437.tif 000179
に変換され、ここで、前記変数ｙは従属変数であり、ｕおよびｒは制御変数であり、
前記反復関数がさらに、前記微分方程式の０≦τ≦Ｔにおける解として効率的な最適化の軌跡を計算するための解法を提供する、請求項１６に記載の最適化システム。
【請求項１８】
ニュートン法を使用してφを前記１次微分方程式
【数９】
2007509437.tif 000180
に変換し、
ここで、偏導関数ベクトルＧは、
【数１０】
2007509437.tif 000181
として定義され、ヘッシアン行列Ｈは、
【数１１】
2007509437.tif 000182
として定義されるサイズＮｘＮの行列である、請求項１７に記載の最適化システム。
【請求項１９】
前記最適化プログラムが、前記最適化プログラムのパフォーマンスを最適化するための制御変数を計算する、請求項１７に記載の最適化システム。
【請求項２０】
前記最適化プログラムが、
前記制御変数を最適化するための最適化問題を定式化すること、
前記最適化変数を計算するための反復を定式化すること、および
第２の連続型反復変数を使用した反復を継続することによって、
前記最適化プログラムのパフォーマンスを最適化するための制御変数を計算する、請求項１９に記載の最適化システム。
【請求項２１】
前記最適化プログラムが、連続最適化問題を常微分方程式に変換することによって前記最適化問題の最小の解に向かう軌跡を求める前記連続最適化問題の解法を採用し、
前記常微分方程式は、連続型反復変数τに関する前記状態ベクトルｙの正規化された変化率を表す以下に示す式であり、
【数１２】
2007509437.tif 000183
ここで、偏導関数ベクトルＧは、
【数１３】
2007509437.tif 000184
として定義され、ヘッシアン行列Ｈは、
【数１４】
2007509437.tif 000185
として定義されるサイズＮｘＮの行列であり、
前記常微分方程式は、代替的に、
【数１５】
2007509437.tif 000186
と表され、ここで、
ｗ(τ)＝Ｈ−１Ｇ(ｙ(τ),ｕ(τ),ｒ(τ))、および
ρ(ｙ(τ),ｕ(τ),ｒ(τ))＝ＧＴＨ−１Ｇ(ｙ(τ),ｕ(τ),ｒ(τ))
である、請求項１３に記載の最適化システム。
【請求項２２】
前記最適化プログラムが、前記最適化問題をｕ(τ)，ｒ(τ)，τ∈［０，Ｔ］に関する制御問題
【数１６】
2007509437.tif 000187
で具体化することによって、準最適方式で軌跡を求める制御変数を計算するための反復法を提供する制御問題に前記最適化問題を変換し、
ここで、Ｑは対称正定値行列であり、前記制御問題は、以下に示す制約条件
【数１７】
2007509437.tif 000188
に従って、前記目的関数Ｆを最適な状態ベクトル変化の軌跡ｙ＊(τ)内の局所的最小値にするような最適制御ポリシーｕ＊(τ)およびｒ＊(τ)を定める、請求項２１に記載の最適化システム。
【請求項２３】
最適な状態ベクトル変化の最適化軌跡ｙ＊(τ)を検出する前記最適制御ポリシーｕ＊(τ)およびｒ＊(τ)（τ∈［０，Ｔ］）が、ハミルトニアン
【数１８】
2007509437.tif 000189
から獲得され、
該ハミルトニアンは、ｕおよびｒの関数としてｕ＊(τ)およびｒ＊(τ)で最小値をとり、
前記状態関数ｙ(τ)および共役関数ｐ(τ)の導関数が、すべてのτ∈［０，Ｔ］について、
【数１９】
2007509437.tif 000190
であり、境界条件、
【数２０】
2007509437.tif 000191
を満たすように、連続型共役関数ｐ(τ)が存在する場合に、ｕ＊(τ)は最適となる、
請求項２２に記載の最適化システム。
【請求項２４】
前記最適化プログラムが、前記最適化問題を、準最適方式で前記軌跡を求める制御変数を計算するための反復法を提供する制御問題に変換し、前記制御問題は、
【数２１】
2007509437.tif 000192
として要約され、ここで、ｙ＊(τ)は準最適状態ベクトル軌跡、ｙ０は初期状態ベクトル、ｐ(τ)は共役関数、ｕ＊(τ)およびｒ＊(τ)は最適制御ポリシー、Ｑは対称正定値行列、τ∈［０，Ｔ］である、請求項２２に記載の最適化システム。
【請求項２５】
前記最適化プログラムは、準最適軌跡を繰り返し計算して前記最適化問題の準最適解を得る、請求項２３に記載の最適化システム。
【請求項２６】
最適化問題の当初の準最適解を使用して、修正および展開された連続型準最適解を生成する修正機能をさらに備える、請求項１に記載の最適化システム。
【請求項２７】
前記修正機能が、
前記当初の準最適解を時間的に展開することと、
最適化問題の前記当初の準最適解が計算されたとき以降に発生したイベントを最適に考慮することと、
特定の時間区間にわたって前記最適化問題の前記当初の準最適解の特定の状態変数を任意選択で固定することと、
前記最適化問題の前記当初の準最適解を計算するために当初使用された決定変数の値と制約を任意選択で変更することと、
最適化問題の前記当初の準最適解に関する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することと、
によって、最適化問題の前記当初の準最適解を使用して、修正および展開された連続型準最適解を生成する、請求項２６に記載の最適化システム。
【請求項２８】
前記当初の準最適解を時間的に展開することが、
時刻ｔ１−ΔＴ＋Ｔからｔ１−ΔＴ＋Ｔまで長さΔＴの時間区間にわたって前記当初の準最適解を時刻ｔ１で展開するために、前記当初の準最適解の前記制御関数を使用し、前記当初の準最適解の前記状態関数の偏微分方程式を積分することをさらに含む、請求項２６に記載の最適化システム。
【請求項２９】
最適化問題の前記当初の準最適解が計算されたとき以降に発生したイベントを最適に考慮することが、
時刻ｔ１−ΔＴから時刻ｔ１までのΔＴの時間区間に発生した前記イベントのすべてを、時刻ｔ１で発生した単一イベントとして近似し、前記単一イベントを使用して、前記修正および展開された連続型準最適解を計算するときに前記当初の準最適解の初期条件を変更することをさらに含む、
請求項２６に記載の最適化システム。
【請求項３０】
特定の時間区間にわたって前記最適化問題の前記当初の準最適解の特定の状態変数を固定することが、
時刻ｔ１から時刻ｔ１＋Ｔまでの特定の時間区間にわたって固定する状態変数を識別することと、
識別された状態変数のそれぞれについて、前記修正および展開された準最適解のための前記識別された状態変数について計算された値が、前記最適化問題の前期当初の準最適解について計算された値と同じになるように、前記修正および展開された連続型準最適解を計算することと、
をさらに含む請求項２６に記載の最適化システム。
【請求項３１】
前記最適化問題の前記当初の準最適解を計算するために当初使用された決定変数の値および制約を変更することが、
フェンス区間にわたって変数を固定することをさらに含む、請求項３０に記載の最適化システム。
【請求項３２】
前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
時間区間Ｔにわたってフェンスのない変数の変動を最小化することをさらに含む、請求項３１に記載の最適化システム。
【請求項３３】
前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
前記最適化問題のノミナル解を時刻ｔ１＋Ｔから時刻ｔ１＋Ｔ＋ΔＴまで分析的に展開することをさらに含む、請求項３０に記載の最適化システム。
【請求項３４】
前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
修正のために特定の区間に関連する変数の優先順位を付けることをさらに含む、請求項３３に記載の最適化システム。
【請求項３５】
前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
増分最適化問題の解を展開する間に基準の制約を変更することをさらに含む、請求項３０に記載の最適化システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、さまざまな問題に対する相対的に最適な解を決定するため、および複雑な問題の解を求める計算プロセスと計算アプローチの構築を制御するために使用する一般最適化問題に関し、特に、一般的な目的関数のクラスを最小化するため、時間の経過につれて解を調整するため、およびさまざまな計算プロセスに対する相対的に最適な計算制御を構築し、実現するための一般的な計算の方法およびシステムに関する。
【背景技術】
【０００２】
本出願は、いずれも２００２年１０月２３日に出願された仮出願第６０／４２０，９２１号および６０／４２０，９２２号の利益を主張する２００３年１０月２３日に出願された米国特許出願第１０／６９３，７２９号の一部継続出願である。
【０００３】
現在、数学および計算による最適化の分野には、膨大な研究および開発の労力が費やされている。最適化の技術は、複雑な系および組織向けの、費用対効果が高く、時間的に効率のよい制御戦略を含めて、さまざまな日常の問題および複雑な理論上の問題に適用され、こうした問題を解決している。たとえば、大規模なメーカー組織における労働力資源と設備資源のスケジューリングに最適化の技術を適用することによって、指定された量の製品を費用対効果が最も高い方法で生産し、さまざまな時点で計画された製品需要を満たすことができる。別の例として、都市の信号機制御のように大規模な制御問題に最適化の技術を適用することによって、都市全体にわたってさまざまな交差点での過度な遅延の発生をできるだけ抑えながら膨大な量の交通を処理できるようにし、さまざまな２次的および３次的な経路を利用して通行できるようにする。多くの自然体系は、物理的な特性と原理によって制御され、多くの物理的な制約によって制限された複雑な定義による問題の相対的な最適解を求めていると考えることができる。たとえば、複雑な化学反応は、分子などの個々の要素が最小のエネルギー状態と最大のエントロピーを追求するという比較的単純な熱力学の原理によって制御されており、分子が量子力学に従ってそれ自体を再構成し、それ自体の間でエネルギーを交換できるさまざまな方法によって制約されると考えることができる。さらに複雑な生体系では、化学と熱力学の原理だけでなく、情報量と情報交換、区分化とモジュール設計、および他に考慮すべき問題によっても制御されると考えることができる。このように、最適化問題には、自然現象や、複雑な系の設計、体系化、および制御や、情報の体系化、送信、および処理を、数学および計算によって表現した極めて多種多様なモデルが含まれる。
【発明の開示】
【発明が解決しようとする課題】
【０００４】
数理的にモデル化された問題の準最適解（near-optimal solution）および最適解を求めるための現在のアプローチは限定されている。単純で典型的な問題で標準的に使用する少数の変数と制約から、現実の世界で一般的に要求される多くの変数と制約まで、こうした問題の決定変数と制約の数が増大すると、準最適解および最適解を求めるために必要な計算のリソースは多くの場合に指数関数的に増大する。現在の技術は、任意の問題に適用するには不十分であり、非常に単純なタイプの問題にしか適用できない。現在使用可能な多くの技術には、変数が連続でなければならないことや、最適化問題の可能性のある解のセットを表す超次元体が必要なこと、問題の領域（domain）が凸（convex）でなければならないことなど、最適化モデルに適用する苛酷な制約が含まれる。
【０００５】
したがって、さまざまなタイプの複雑な系の挙動に関する研究者、開発者、システム・モデル作成者、および調査員は、より一般的に適用できる効率的な最適化の方法およびシステムの必要性を認識してきた。
【課題を解決するための手段】
【０００６】
本発明のさまざまな実施形態には、準最適解または最適解を求める問題を数理的にモデル化する高次元領域の関数を最小化するアプローチによって問題を最適化する一般的な方法およびシステムが含まれる。このような実施形態は、記号形式で（symbolically）システムの数理的な記述を受け取る。こうした記述には、実数値、整数値、ブール値（Boolean value）など、さまざまな型の決定変数（decision variables）が含まれる。準最適解または最適解を求めるために最小化される高次元領域の目的関数は、不等式制約（inequality constraint）と等式制約（equality constraint）を含む、決定変数の値に関するさまざまな制約を伴う。本発明のさまざまな実施形態では、目的関数および制約が広域的な目的関数に組み込まれている。この関数の局所的最小値（local minima、極小）および広域的最小値（global minima）は、臨界点（critical points）で発生する。こうした臨界点の１つまたは複数は、モデル化された問題の準最適解または最適解である。本発明のさまざまな実施形態では、広域的な目的関数と臨界点を検出する手続きとが、連続変数およびパラメータによる微分方程式系に変形されるので、微分方程式を解くための強力な多項式時間法（polynomial-time method）を適用して関数の臨界点を特定できる。これらの実施形態は、内点に基づく方法（interior-point-based method）を実行しており、広域的な勾配降下フィールドの定式化を使用して広域的な目的関数の超次元領域内を効率的に移動し、局所的最小値を検出する。局所的な臨界点の近くでは、広域的な降下フィールドは変数のより良い近似を提供する場合に効果的ではないが、本発明の実施形態が提供する操作方法では、局所的な臨界点から解に至る局所的な降下フィールドを生成する。本発明のさまざまな実施形態では、局所的な臨界点の近傍を過ぎると、広域的な勾配降下法が再開され、超次元領域を経由して局所的最小値に至る降下のパスを再びたどるための適切な方向を確定する方法を使用する。こうした局所的最小値の１つは、問題の準最適解を示している。本発明の実施形態では、広域的な勾配降下フィールドの最適化法と複数のエージェントを使用した局所的な勾配降下フィールドの最適化法とによる分散および分解も可能であり、パラレルな計算と時間的な効率向上を実現している。本発明のさまざまな実施形態には、最適化問題の解の時間に関してかなり連続的に調整するためのアプローチがさらに含まれており、さまざまなイベント、優先順位の変化、および予測の変化に対応するので、新しい最適解を再計算する必要はない。本発明の多くの実施形態は、特にさまざまなクラスの最適化問題を対象とするが、本発明の他の実施形態は複雑な階層型の計算プロセスを構築するため、および一般的な計算プロセスを最適またはほぼ最適に制御するためのより一般的なアプローチを提供する。
【発明を実施するための最良の形態】
【０００７】
本発明は、数理的に表現された問題の一般的な最適化、最適制御問題の解の比較的連続する調整、およびより抽象的には計算プロセスを最適またはほぼ最適に制御するための高レベルの階層構造に関する。本発明のさまざまな実施形態には、単一プロセッサのコンピュータ・システム上で動作する、またはマルチプロセッサ・コンピュータ・システム上で並列で、または多くの分散型の相互接続された単一プロセッサおよび／またはマルチプロセッサ・コンピュータ・システム上で動作する複雑なソフトウェア・プログラムが含まれる。このようなプログラムは、完全に数理的な特性を備えており、複雑で記述の困難な数理的技術を必要とする。こうした理由により、本発明について、具体的な問題に関連付け、多くの数式を使用して、さらに多くの図面と制御流れ図を参照しながら、以下でその一部について説明する。本発明のさまざまな実施形態を表すプロセスを実施するプログラムは、ソフトウェア開発技術者に対しては数式のみで十分に説明でき、特徴付けできるが、以下の説明で示すより図式的で具体的な問題指向の制御流れ図によるアプローチは、本発明のさまざまな実施形態を個別の経歴をもつ読者に特にわかりやすいように、さまざまな方法で説明することを目的としている。本発明について、以下の項、（１）最適化、（２）修正、（３）メタ制御とプロセス階層構造で説明する。
【０００８】
最適化
図１Ａ〜Ｂは、非常に単純な２次元の最適化問題を示している。図１Ａと１Ｂは２次元面の２つの異なるビューを示しており、この面は数理的に次のように表される。
ｚ＝ｆ(ｘ１，ｘ２)
ここで、問題の領域は２次元面であり、問題の目標は従属変数ｚが最小となるｘ１およびｘ２の値を求めることである。換言すると、面のグラフは２つの変数、ｘ１およびｘ２の関数によって特徴付けられ点（ｘｌ，ｘ２，ｚ）の集合を、ｘ１軸１０２、ｘ２軸１０４およびｚ軸１０６についてプロットしたものである。具体的な例として、ｚは長さｘ１および幅ｘ２の装置の製造に関連するコストを表すことができる。このモデルに関する最適化問題は、装置を製造するコストを削減し、可能な限り最小の値とする最適な装置の長さと装置の幅を検出することでもよい。多くの場合には、最適化問題の解を求めるための閉形式の表現も、数理モデルさえも利用できない。代わりに、決定変数のいくつかの実現可能な初期値（この場合は２次元面１１０上の点に対応するｘ１およびｘ１の値）で開始し、実現可能な初期値１１０から関数が局所的最小値値または広域的最小値値ｚｍｉｎをとる座標（ｘ１ｍｉｎ，ｘ２ｍｉｎ，ｚｍｉｎ）である面上の点１１４に至るパス（図１Ａに濃い曲線１１２で示されている）を検出しようとする。
【０００９】
ただし、問題に対する準最適解または最適解を求めるためのアプローチの１つは、実現可能なすべての（ｘ１，ｘ２）のペアについてｚを計算するものであることに留意されたい。しかし、一般にこのようなアプローチは非常に単純な問題でも計算コストがかかる。多くの場合に、決定変数ｘ１およびｘ２などの決定変数の領域は実数の部分集合または他の自然領域の部分集合である。このことによって完全な組合せ検索を行うために必要な関数を評価する決定変数の実現可能な値の数は大幅に減少する。しかし、このような場合でも、組合せ検索法は一般的に実現可能ではない。
【００１０】
図１Ａと１Ｂのような２次元面の図を以降の説明で使用するが、これらは説明のみを目的として示されていることに留意されたい。現実の世界の問題では、数千、数十万、数百万の変数と制約が存在するため、極めて高次元の問題領域になる場合がある。一般に、このような問題領域は、制約によって指定される超次元面で囲まれた超次元体または超次元多様体とみなすことができる。このような問題領域を図示することは不可能であるが、超次元の問題領域に対処するための手法は３次元の図から類推できる。たとえば、２次元の問題領域における実現可能点は２次元面の点であるのに対して、超次元の問題領域における実現可能点は超次元境界面で囲まれた超次元体または超次元多様体の内部の点である。
【００１１】
数理モデルの領域内で、最初の実現可能点から局所的最小値点または広域的最小値点まで移動するために、その他の多くのアプローチを利用できる。たとえば、１つのアプローチでは、各点で関数の勾配を計算し、面に沿って最も急勾配で降下する方向に進み、局所的最小値または広域的最小値に至る。このアプローチは、唯一の広域的最小値が存在する凸面の問題には最適である。実世界の最適化問題の超次元体領域は、凸面上でない場合が非常に多く、広域的最小値以外に、鞍点、極大、局所的最小値を含む多くの臨界点を含んでいる。したがって、急降下の勾配を使用したアプローチは超次元の問題領域で局所的最小値に移動する場合には使用できるが、複数の局所的最小値を操作し、評価するためには検索の手法を使用しなければならない。たとえば、図１Ａでは、最急降下法によって決定された軌跡は点１１６で示される局所的最小値で終わっており、点１１２の広域的最小値を検出し、評価するには他のいくつかのアプローチが必要である。図１Ｂは問題領域の面と広域的最小値へのパスを示しており、図１Ａで示されたものをｚ軸に関して約４５度回転させたものである。
【００１２】
図２は、本発明の１つの実施形態を表す最適化の方法およびシステムの最高レベルのビューを示す制御流れ図である。最初のステップ２０１で、最適化される特定の系または問題が数理的にモデル化され、その数理モデルが後続の最適化のステップへの入力となる。本発明の多くの実施形態には、問題の数理モデルの作成を容易にするためのグラフィカル・ユーザ・インタフェース、テキストベースのインタープリタ、および他のインタフェースが含まれる。以下に、仮説的な問題の特定の数理モデルを、最適化問題の数理モデルの例として示す。一般に、数理モデルは記号形式（symbolic form）で提供され、シンボリックな項と演算子を含む数式で構成される。次に、ステップ２０２で、本方法およびシステムは前処理のステップを実行し、記号形式で表現された数理モデルから本発明の実施形態を表す方法で必要とする標準形式（standard form）に変換する。前処理のステップについては以下で仮説的な数理モデルに関連して説明するので、ここでは前処理のステップについて形式的に説明する。次に、ステップ２０３では、数理モデルの記号形式による標準の形式が使用され、繰り返しを設定し、内点法の例と勾配降下フィールド（descent field）によって導かれた軌跡を使用して得られた関数を繰り返しによって最適化し、超次元問題領域内にある１つ以上の局所的最小値を識別する。降下フィールド・メカニズムのさまざまな戦略を使用して、存在する可能性のある多数の局所的最小値を検出して評価し、数理的にモデル化された問題の準最適解を特定する。本発明のさまざまな実施形態では、さまざまな検索戦略を採用している。本発明の実施形態を表す方法でも、現在採用されている他の最適化法でも、指定された許容範囲内または妥当な時間内で、妥当な量のコンピュータ・リソースを使用して、広域的最小値を必ず検出することを保証しているわけではないことに留意するのは重要である。しかし、多くの実際的な最適化問題には境界が指定されており、局所的最小値の数は限られている。最終的に、指定された量のコンピュータ・リソースを使用して、指定された許容範囲または最大時間内で最善の解が見つかると、後処理のステップ２０４でこの解を変換し、ステップ２０１で初めに提供された問題に対する数理モデルのコンテクストで１つまたは複数の準最適解を提供する。
【００１３】
数理的なモデル化のステップと前処理のステップ、すなわち図２の２０１と２０２は、それぞれ具体的、仮説的な問題を使用しておそらく最もよく理解される。図３は、数理的なモデル化のステップと前処理のステップを説明するために使用する仮説的な問題を示している。この仮説的な問題には、ｐ個の製造工場３０２〜３０４、ｗ個の卸売店３０６〜３１１、ｃ個の顧客３１４〜３１８が存在する。ある製造工場、たとえば製造工場３０２は、図３のさまざまなタイプの線（たとえば製造工場３０２と卸売店３０６を結ぶ破線３２０）で示されるように、ｗ個の卸売店のいずれにも商品を出荷できる。さらに、各卸売店は、さまざまなタイプの線、たとえば卸売店３０６と顧客３１４を結ぶ破線３２２で示される卸売店と顧客の間のパスに従って、ｃ個の顧客のいずれにも商品を出荷できる。図３には、さまざまなスタイルの線を使用して、特定の製造工場からさまざまな卸売店へのパス、および特定の卸売店からさまざまな顧客へのパスが示されている。
【００１４】
図４〜７は、仮説的な問題のその他の側面と特性を、図３の描画法を使用して示している。たとえば、ｘＰＷｉｊの表記は製造工場ｉから卸売店ｊに出荷された商品の量を示し、ｘＷＣｊｋの形の項は卸売店ｊから顧客ｋに出荷された商品の量を示している。ｗ個の卸売店の候補が存在するが、ｗ個すべての卸売店を使用する必要はない。ブール変数ｂｊは卸売店ｊが使用されているかいないかを示し、ブール値「１」すなわち「ＴＲＵＥ」は卸売店が使用されていることを、ブール値「０」すなわち「ＦＡＬＳＥ」は卸売店が使用されていないことを示す。図４では、たとえば、使用されている卸売店を示す４角形の実線の枠と使用されている卸売店に関連付けられられたブール値「１」で示されるように、卸売店１（３０６）、３（３０８）および５（３１０）が使用されており、卸売店２（３０７）、４（３０９）、および卸売店ｗ（３１１）を含めて卸売店５より後のすべての卸売店は使用されていない。図４では、使用されている卸売店は括弧内に連番で別の番号が付けられており、３つの使用されている卸売店には、順に１、２、およびｎの番号が付けられている。ただし、ｎは使用されている卸売店の数であり、図４の例では３である。
【００１５】
図５は、製造工場から卸売店に商品を出荷する際のコスト、および卸売店から顧客に商品を出荷する際のコストを示している。ＣＰＷｉｊの形の項は製造工場ｉから卸売店ｊに商品を出荷する品目単位のコストを示し、ＣＷＣｊｋの形の項は卸売店ｊから顧客ｋに商品を出荷する品目単位のコストを示している。図５では、出荷ルートのそれぞれについて、対応する品目単位のコストが表示されている。
【００１６】
図６は、仮説的な問題の付加的な特性を示している。各製造工場には製造能力が関連付けられており、製造工場ｉの出荷能力を示すｓｉの形の項がある。各卸売店には受け入れ能力があり、これは卸売店ｊが商品の受け入れ能力を示す項ｖｊで表される。また、各卸売店には運用コストがあり、卸売店ｊの運用コストを示すｒｊの形の項がある。最後に、各顧客には商品の需要が関連付けられている。ｄｋの形の項は、顧客ｋの商品に対する需要を示す。図６には、製造工場、卸売店、および顧客に関連する供給能力、受け入れ能力、運用コストおよび需要の詳細がすべて明確に示されている。最後に、図７は２つの製造工場、３つの潜在的な卸売店、５つの顧客からなる単純な仮説的問題を示しており、関連のコスト、能力および需要を数字で示している。
【００１７】
図３〜７に示す最適化問題は、顧客のすべての需要を満たし、能力が超えない範囲での卸売店の運用、製造工場から卸売店への出荷量、卸売店から顧客への出荷量に関する最小限のコスト構成を算出することである。最適化問題のより形式的な記号形式による数理的な表現は以下のとおりである。
【００１８】
まず、特定のシステム構成に関連する総コストは次の式で表される。  
【数２２】
2007509437.tif 000002
つまり、総コストは製造工場から卸売店、卸売店から顧客へのすべて出荷に関連するコストと、卸売店の運用コストの合計である。この仮説的な最適化問題において、総コストの最小値を求める。さらに、このモデルにはさまざまな制約が含まれる。まず、製造工場は製造工場の能力を超える量の商品を出荷することはできない。これは次のように表される。  
【数２３】
2007509437.tif 000003
【００１９】
すべての顧客は、少なくともその需要に等しい量の商品を受け取る。これは次のように表される。  
【数２４】
2007509437.tif 000004
【００２０】
さらに、製造工場から各卸売店に出荷する商品の量は、卸売店の商品受け入れ能力以下である。これは次のように表される。  
【数２５】
2007509437.tif 000005
同様に、卸売店から出荷する商品の量は、卸売店の商品受け入れ能力以下である。これは次のように表される。  
【数２６】
2007509437.tif 000006
卸売店には商品を保管しない。各卸売店は、受け取った商品をすべて出荷する。これは、次に示す保管に関する制約によって表される。  
【数２７】
2007509437.tif 000007
【００２１】
図３〜７に示す出荷のパスは単方向である。つまり、商品は返品されない。これは次のように表される。  
ｘＰＷｉｊ≧０
ｘＷＣｊｋ≧０
仮説的な問題では、顧客の需要を満たすのに十分な製造能力と卸売店の受け入れ能力があることが仮定される。これは次のように表される。  
【数２８】
2007509437.tif 000008
商品は、運用している卸売店からのみ出荷できる。これは次の２つの制約によって表される。  
ｘＰＷｉｊ(１−ｂｊ)＝０
ｘＷＣｊｋ(１−ｂｊ)＝０
運用している卸売店の数は、卸売店の候補の総数ｗ以下である。これは次のように表される。
【数２９】
2007509437.tif 000009
最終的に、項ｙｊはブール値（Boolean value）ｂ∈{０，１}である。
【００２２】
前述の数理モデルは、コストの関数と関連の制約をより簡潔に表すように若干変形してもよい。まず、出荷する商品の量を表すｘＰＷｉｊおよびｘＷＣｊｋのタイプのパラメータを、次のようにスケーリングされた出荷の項で置き換えてもよい。  
ｍＰＷｉｊ＝ｘＰＷｉｊ／ｓｉ、０≦ｍＰＷｉｊ≦１
ｍＷＣｊｋ＝ｘＷＣｊｋ／ｓｋ、０≦ｍＷＣｊｋ≦１
前述の卸売店の能力に基づく制約は、冗長になるので制約のセットから削除でき、問題のモデルまたはモデルに関して得られる解に影響することはない。卸売店に商品を保管しないことを表す保管の制約は、スケーリングされた出荷の項を使用するために次のように変更される。  
【数３０】
2007509437.tif 000010
運用されていない卸売店からは商品を出荷しない。これは次の制約によって表される。  
ｍＷＣｊｋ≦ｂｊ
総コストは、商品のスケーリングされた出荷の項を使用して次のように表すことができる。  
【数３１】
2007509437.tif 000011
能力の制約は、商品のスケーリングされた出荷の項を使用して次のように表される。  
【数３２】
2007509437.tif 000012
前述の能力に関する制約は、次のように、０に関する不等式制約として作成し直すことができる。  
【数３３】
2007509437.tif 000013
ここで、運用している卸売店の数は、卸売店の候補の総数以下でなければならない、商品は返品できないこと、さらにｂｊ項はブール値であるという事実は次のように表される。
【数３４】
2007509437.tif 000014
【００２３】
図８は、図２に示す高レベルの制御流れ図のステップ２０２で呼び出される前処理のルーチンを示すより詳細な制御流れ図である。図８に示す前処理ルーチンの各ステップについては、前述の仮説的な数理モデルの対応する変換に関連付けながら以下で詳細に説明する。まず、ステップ８０１で、前処理ルーチンは、制約を含めて最適化問題を最初に数理的に表現した記号形式のモデルを受け取る。ステップ８０３で、前処理ルーチンは冗長な制約をできるだけ除去するように試みる。冗長な制約の除去の例は、上の変換されたモデルで実行されている。つまり、元の形のモデルにおいて冗長であると認識された卸売店の能力に関する制約は、変換されたモデルからは削除されている。次にステップ８０５で、すべてのブール値および整数値の決定変数（たとえば上の仮説的な例におけるブール変数ｂｊ）は、浮動小数点変数に変換される。ステップ８０７で、結果として得られる数理的に表現された記号形式のモデルは標準の形式に配列される。ステップ８０９で、等式制約はゆるやかな不等式制約に置き換えられる。各等式制約は、正の値をとる追加の変数を使用する３つの不等式制約に緩和される。ステップ８１１で、追加の変数ｚと追加の制約が追加され、基準が追加の制約に変換される。最後にステップ８１３で、項を共に加えることによって、変換された最適化問題を内点法（Interior Point Method）のバリア関数（Barrier function）に変換する。各項は、制御変数をかけた制約を表す式に関するバリア関数の評価である。内点法のバリア関数には、２つの制御変数のセット、ｕとｒが含まれる。前処理ルーチンの結果として、数理モデルは、図２のステップ２０３に示す反復最適化の技術によって最適化できる形の記号形式で表現される。
【００２４】
図９は、図７に示す特定の最適化問題を最初に数理的に表現した記号形式のモデルを示している。コスト関数９０２は、整数変数ｘおよびブール変数ｂのベクトルに関して最適化される。制約には、製造工場の能力に関する制約９０４、顧客の需要に関する制約９０６、保管に関する制約９０８、最適化で選択されていない卸売店には出荷も入荷も発生しない、すなわち対応するブール変数はゼロであるという事実をモデル化した制約９１０、３つの卸売店のうち２つだけが運用されることを示す制約９１２、顧客から卸売店へも卸売店から製造工場へも商品の積み戻しがないことを示す制約９１４、およびｂｊがブール値をとる変数であることを示す制約９１６が含まれる。
【００２５】
前処理ルーチンの第３のステップである図８のステップ８０５には、ブール変数、整数変数、およびその他の浮動小数点以外の変数を浮動小数点変数に変換するステップが含まれる。各ブール変数ｂｉは、範囲［０，１］にあり、追加の制約を伴う浮動小数点の値ｙｉに置き換えられる。  
ｙｉ（ｌ−ｙｉ）＝０
追加の制約は、浮動小数点変数ｙｉが０または１の２つの値のうち、いずれか１つしか取ることができないことを示している。各整数変数は、ｋ個の項からなる等価な２進多項式（base-2 polynomial）に置き換えられる。ここでｋは整数変数がとり得る最大値の２を底とする対数にほぼ等しい。２進多項式の各項は２のべき乗にブール値の係数をかけたもので構成される。したがって、２進多項式は［０，１］のすべての範囲で次に示す形のｋ個の制約を伴うｋ個の浮動小数点変数ｙ０，ｙ１，．．．ｙｋ−１で置き換えられる。  
ｙｉ（ｌ−ｙｉ）＝０
さまざまな型の決定変数の浮動小数点数への変換は、形式的に次のように要約される。  
【数３５】
2007509437.tif 000015
【００２６】
図１０は、仮説的な数理モデルから浮動小数点変数のみを使用したモデルへの変換を示している。
【００２７】
次に、数理モデルが標準の形式に配列される。標準の形式によるシンボリックな問題の表現は、形式的には次のように定義される。  
【数３６】
2007509437.tif 000016
つまり、標準の形式のモデルには、コストの関数Ｃ（ｘ）、Ｍ１個の不等式制約
【数３７】
2007509437.tif 000017
および複数個の等式制約  
【数３８】
2007509437.tif 000018
が含まれる。  
図１１は標準の形式に変換した仮説的な問題の数理モデルを示している。元の等式制約と不等式制約を標準の形式で必要な右辺が０となる制約に変換するのは、まだ標準の形式になっていない制約式の両辺に項を加えるだけの簡単な作業であることに留意されたい。
【００２８】
前処理ルーチンの次のステップ（図８の８０９）は、すべての等式制約を不等式制約に変える操作である。これを行うために、各等式制約を２つの不等式制約に変換する。各不等式制約には、元の等式制約の左辺の式と正の束縛変数（binding variable）が含まれる。等式制約ｈｊ（ｘ）から２つの不等式制約への変換は、形式的に次のように表される。  
【数３９】
2007509437.tif 000019
さらに、制約
ｓｊ≧０
が追加され、ｈｊ（ｘ）を表す制約が実現可能であることを保証する。
【００２９】
したがって、以下の等式制約から不等式制約への変換の後で、制約の総数は問題の標準の形式では元の等式制約の数に対して２ｋだけ増える。図１２Ａ〜Ｄは、変換された等式制約における束縛変数の効果を表している。図１２Ａは、値０を持つ元の等式制約ｈｊ（ｘ）を示している（１２０２）。前述のように等式制約ｈｊ（ｘ）が３つの不等式制約に変換されると、ｈｊ（ｘ）の値の範囲は０から−ｓｊ  １２０４とｓｊ  １２０６の間にあるすべての値に拡張される。したがって、等式制約を３つの不等式制約に変換することによって、実質的に等式制約は点から区間に緩和される。後述するように、等式制約を緩和すると、等式制約によって１つまたは複数の局所的最小値に向かう最初の降下を過剰に拘束したり支配したりすることはなくなるので、反復最適化法が円滑化される。また、等式制約によって最適化の軌跡は問題領域の境界になることを余儀なくされるので、等式制約を緩和すると内点降下も促進される。しかし、最適化が収束し始めたら、等式制約の緩和を排除する必要がある。これは図１２Ｃおよび１２Ｄに示すように、バリア変数ｓｊの絶対値を減少し、ｈｊ（ｘ）の値に関する制約をさらに強めて０を中心にして対称的に配置される小さな領域に拘束することで達成できる。等式制約の緩和の程度は、基準に加えられた各束縛変数ｓｊの係数ｒｉで表すことができる。等式制約の不等式制約への変換は、形式的に次のように表される。  
【数４０】
2007509437.tif 000020
【００３０】
図１３は、仮説的な問題の等式制約を不等式制約に変換した後の数理的に表現された記号形式のモデルを示している。図１４は、束縛変数ｒｉの効果を図形的に示している。図１４では、最適化法は初期点１４０２から始まり、カップ型の面１４０６で囲まれた問題領域１４０４内を最小１４０８に向かって降下する。最初は、束縛変数ｓｉの係数ｒｉは小さく、図１４において円盤形の領域１４１０で示されるように、初期点１４０２に関する値の範囲全体は大きい。問題領域を降下して最小１４０８に近づくにつれて、最適１４０８の近くでは、束縛変数ｓｉの係数ｒｉは大きくなり、解の集合に対する制約はしだいに厳格になる。最適化法が進むに従ってｓｉはどんどん小さくなるので、変換された等式制約は再び極端に厳しく制約される（１４１２）。
【００３１】
前処理ルーチンの次のステップには、問題の数理モデルにエンベロープ変数（envelope variable）ｚと追加の制約を追加する処理が含まれる（図８の８１１）。エンベロープ変数ｚは、次のように定義される。  
【数４１】
2007509437.tif 000021
したがって、次の不等式制約が得られる。  
【数４２】
2007509437.tif 000022
エンベロープ変数ｚを導入した後、問題の形式的な記述は次のようになる。  
【数４３】
2007509437.tif 000023
【００３２】
図１５は、エンベロープ変数ｚを導入した後の仮説的な問題の数理モデルを示している。
【００３３】
図１６Ａ〜Ｄはバリア関数の概念を示している。たとえば、最適化問題の問題領域が図１６Ａに示すように放物線に基づくくさび形１６０２であると仮定する。最適化の初期点１６０４では相対的に大きなｚ値１６０６をとり、最小１６０８で０に減少する。図１６Ｂは、問題領域１６０２のｘ−ｙ平面への投影図１６１０を示している。初期点１６０４から最小点１６０８までの最適化の軌跡は３つの線分１６１２〜１６１４で示されている。内点法では、問題領域を移動するとき、最適化の軌跡は問題領域内にとどまる必要がある。換言すると、問題領域内で最小に至る降下の軌跡上にある現在の点は、常に問題領域内すなわち超次元問題領域内の多様体の境界の内部に残っている必要がある。たとえば、図１６Ｂで示されるように、最初の軌跡１６１２は問題領域に対応する形状の境界に向かう角度で進み、点１６１６で境界と交わる。引き続き同じ方向に進むと、境界を越えて破線の線分１６１８に沿って問題領域外のある点を通過することになる。最適化の軌跡を問題領域内に維持するには、問題領域内では影響が小さいか全くないが、最適化の軌跡が問題領域の境界面の近くになると、軌跡を境界から問題領域に戻すために値が増大するようなバリア関数を問題のモデルに含める必要がある。たとえば図１６Ｂでは、最適化の軌跡１６１２が点１６１６に近づくに従って、数理モデルに追加されたバリア関数が急速に増大し始め、エンベロープ変数ｚの値を増大することによって最適化の軌跡を問題領域の内側の境界面から離れるように屈曲させる。屈曲した後は、後続の最適化によって境界点１６２０に向かう新しい第２の軌跡に続く最適化の軌跡１６１３が確立される。境界点１６２０の直前でバリア制約によって再び目的関数ｚの値が大きくなって再び最適化の軌跡の向きが問題領域の内側に屈曲し、新しい最後の軌跡１６１４を得る。
【００３４】
図１６Ｃ〜Ｄはバリア関数の動作を示している。図１６Ｃに示すように現在の最適化の軌跡１６２２を想定し、図１６Ｄに示すように、図１６Ｃに示す現在の軌跡１６２２に沿って問題の領域内を移動する目的関数ｚのｘ軸に関する値について考察する。図１６Ｄからわかるように、ｚの値は現在の軌跡１６２２に沿った問題領域の内部では比較的小さい。しかし、ｘが減少して左の境界上の点１６２４に近づくか、増大して右の境界上の点１６２６に近づくと、バリア関数が影響してそれぞれのｚの値１６２８と１６３０が急速に増大する。換言すると、軌跡を問題領域の内部に拘束するために、モデルに含まれる不等式のそれぞれについて、内点の目的関数にバリア関数を含む項が追加される。軌跡が問題領域から出た場合、つまり問題領域を囲む面の外側に達した場合は、本発明のいくつかの実施形態を表す最適化法は失敗する。
【００３５】
モデルにバリア関数に含む項を追加するのは、最適化の軌跡が問題領域内に確実にとどまるようにするための簡単な方法である。問題領域を囲む面を指定するのは制約なので、制約ごとにバリア関数が追加される。本発明のさまざまな実施形態で使用する最適化技術に適したバリア関数項は、こうした項に関連付けられた制御変数である正の変数の係数ｕｉに制約の式ｇｉの対数をかけ、負号を付けたものである。つまり、次のように表される。  
−ｕｉlog(ｇｉ)
制約の式ｇｉの値が０に向かって減少するに従って、制約の対数log(ｇｉ)の値は急速に減少し、これに対応して絶対値は急速に増大する。Log(ｇｉ)の負号を付けることで、制約式ｇｉが０に非常に近くなると、急速に増大する値が得られ、極限では無限大に近づく。したがって、制約がさらに減少して０になるのを防止できる。この変数の係数ｕｉは、前述の束縛変数の係数ｒｉと同様に、制御変数とみなすことができる。最適化が局所的最小値に向かって最も効率のよい軌跡をたどるように、最適化の間に係数ｕｉとｒｉが変更される。しかし、係数ｕｉは係数ｒｉと違って初期値が大きく、本発明の１つの実施形態によって最適化の軌跡が局所的最小値に近づくに従って減少する。これは、局所的最小値はしばしば問題領域の境界面に近いか境界面上にあるためである。目的関数にバリア関数を追加する形式的な式は次のようになる。  
【数４４】
2007509437.tif 000024
【００３６】
図１７は、バリアの制約を追加した後の仮説的な問題の記号形式のモデルを示している。バリア関数を追加することによって、最初に説明した実施形態の前処理ルーチンを具体化した前処理のステップは完了する。
【００３７】
次に、本発明の１つの実施形態を表す中心的な最適化法について概説する。図１８は、問題領域内の広域的勾配降下フィールドを使用して、初期点から局所的最小値に至る最適化の軌跡の方向を示すようすを表している。問題領域内の各点、たとえば点１８０２などで、目的関数Ｆ(ｘ，ｓ，ｚ，ｕ，ｒ)の偏導関数ベクトルが、目的関数の決定変数ｘｉ、束縛変数ｓｉおよびエンベロープ変数ｚに関する１次偏導関数のベクトルとして計算される。同様に、決定変数ｘｉ、束縛変数ｓｉ、およびエンベロープ変数ｚに関するＦ（ｘ，ｓ，ｚ，ｕ，ｒ）の２次偏導関数の行列として与えられる、関数Ｆ（ｘ，ｓ，ｚ，ｕ，ｒ）のヘッシアン（Hessian）行列は、領域のすべての点で計算可能である。目的関数の勾配は、領域内の各点においてヘッシアン行列の逆行列とＦ（ｘ，ｓ，ｚ，ｕ，ｒ）の偏導関数の積で得られるベクトルとして定義される。これは、決定変数と束縛変数に関して目的関数の値が降下する方向を表しており、点１８０２からの方向１８０４を示している。勾配は問題領域内の各点で計算できるので、目的関数の勾配は問題領域内における降下方向を示し、最小へのパスを指定するベクトル場を表している。一般に、図１８の点１８０２のように局所的最小値から離れた点の勾配は絶対値が相対的に大きく、点１８０６のように局所的最小値に近い点の勾配は絶対値が相対的に小さくなる。
【００３８】
最適化法の数理的記述は、多くの表記法によって簡便化される。まず、目的関数Ｆ（ｘ，ｓ，ｚ，ｕ，ｒ）は、Ｆ（ｙ，ｕ，ｒ）として表現できる。ここで、
【数４５】
2007509437.tif 000025
である。  
【００３９】
したがって、目的関数は、決定変数、束縛変数、エンベロープ変数からなるベクトルｙ、バリア関数の係数からなるベクトルｕ、および束縛変数の係数からなるベクトルｒの関数と考えることができる。ベクトルｙは「状態ベクトル（state vector）」と呼ばれ、ベクトルｕとｒは「制御ベクトル（control vectors）」と呼ばれる。ベクトルｙはＮ個の要素を持っている。ただし、Ｎ＝ｎ（決定変数の数）＋ｋ（束縛変数の数）＋１（エンベロープ変数の数）である。
【００４０】
図１９は、最適化の軌跡を示している。軌跡１９０２は、ベクトルｒとｕの長さを横軸、目的関数Ｆの値を縦軸としてプロットされている。最適化の軌跡は、初期点１９０４から開始され、局所的最小値に向かって目的関数Ｆの値が降下する。降下する間に、ベクトルｕの長さで表されるバリア関数の値は０に向かって減少し、ベクトルｒの長さで表される束縛変数の係数の値は増加する。軌跡は長さΔτの線分に分割される（たとえば最初の線分１９０６）。最適化の間に、ベクトルｙとｕは目的関数の値を小さくするように継続的に再計算される。図１９では、ベクトルｒで表される束縛変数の係数の値も、各区間Δτの間に再計算されることが示されている。しかし、ここで説明する本発明の実施形態では、束縛変数の係数は最適化の軌跡の終わり近くでのみ調整される。
【００４１】
離散的な最適化の技術は次のように表すことができる。
ｙｋ＋ｌ＝ｙｋ＋φｋ（ｙｋ，ｒｋ，ｕｋ）
ここで、ｋ＝０，１，．．．，ｋconvergenaeである。換言すると、計算による最適化法の繰り返しごとに、現在の状態ベクトルの値ｙｋと、ｙｋ、ｒｋ、およびｕｋに依存する離散関数φｋとから新しい状態ベクトルの値ｙｋ＋１が計算される。ただし、反復回数ｋは０からｋconvergenceまでの値をとり、ｋconvergenceで目的関数が局所的最小値となる点の状態ベクトルを決定できる。最適化の最初のステップでは、状態ベクトルと制御ベクトルｙ、ｕ、およびｒは連続型反復パラメータ（continuous iteration parameter）τを媒介変数として使用する。ここでτの範囲は０から反復限界（iteration horizon）Ｔまでである。連続型反復変数τをベクトルｙ、ｕ、およびｒの媒介変数として使用することにより、前述の繰り返しはｙ変数が従属変数であり、ｕとｒは制御変数である１次元微分方程式に変換される。微分方程式に基づく数理的な周知の技術を使用すると、この微分方程式の解として効率的な最適化の軌跡を計算できる。このことについては、以下で説明する。これで、目的関数は形式的に次のように表される。  
Ｆ(ｙ(τ)，ｕ(τ)，ｒ(τ))
ここで、０≦τ≦Ｔである。たとえば、ｙｋ＋１＝ｙｋ＋φｋ(ｙｋ，ｒｋ，ｕｋ)がニュートン法による反復（Newton iteration）の場合は、対応する微分方程式は次の形である。  
【数４６】
2007509437.tif 000026
【００４２】
ニュートン法（Newton’s mothod）は、ｄｙ／ｄτの式を提供する微分方程式を構築するための方法の１つにすぎないことに留意されたい。この方法は、計算を単純にし、容易にするために選択されている。ただし、本発明の実施形態では、ｄｙ／ｄτを近似するために、他の多くのアプローチのいずれを使用してもよい。
【００４３】
以下の説明では、目的関数Ｆの偏導関数ベクトルＧとヘッシアン行列Ｈを頻繁に使用する。偏導関数ベクトルＧは、次のように定義される。  
【数４７】
2007509437.tif 000027
ヘッシアン行列ＨはサイズＮｘＮの行列であり、次のように定義される。  
【数４８】
2007509437.tif 000028
連続する反復変数τに関する状態ベクトルｙの正規化された変化の率は、ヘッシアン行列ＨとベクトルＧを使用して次のように表現される。  
【数４９】
2007509437.tif 000029
または、次のようにも表される。  
【数５０】
2007509437.tif 000030
ここで、ｗ(τ)＝Ｈ−１Ｇ(ｙ(τ), ｕ(τ), ｒ(τ))、および
ρ(ｙ(τ), ｕ(τ), ｒ(τ))＝ＧＴＨ−１Ｇ(ｙ(τ), ｕ(τ), ｒ(τ))
である。
【００４４】
最適化問題は、一般的な最適制御問題の特殊な事例である最短時間問題に書き換えることができる。換言すると、目的関数Ｆ（ｙ，ｕ，ｒ）を連続にすることによって、連続型方法を使用して、目的関数Ｆ（ｙ，ｕ，ｒ）の値を初期点から最小に導く状態変数の変化で構成される最適化の軌跡ｙ（τ）を計算するための計算のアプローチが得られる。しかし、問題領域が非常に大きく、超次元であるために、詳細な計算上の制約またはメタレベルの制約がないと、問題領域内の降下勾配フィールドの方法が得られても、計算による勾配フィールドの降下には収束が遅いという問題がある。
【００４５】
最適制御問題を定式化する。このことによって、収束限界を最小化し、制御変数に関する要件を満足することによって、収束を速めようとする。目的関数Ｆ（ｙ（τ），ｕ（τ），ｒ（τ））の最小化を制御するために、最適制御の定式化は、制御変数ｕ（τ），ｒ（τ）を使用し、収束限界を最小化する制御変数の値を選択することによって、最適化の軌跡を求める。最適化の軌跡を制御する最適制御問題は、ｕ（τ），ｒ（τ），τ∈［０，Ｔ］に関して形式的に次のように表現される。  
【数５１】
2007509437.tif 000031
ここで、Ｑは対称正定値行列（symmetric positive definite matrix）であり、最適制御問題は、次の制約条件に従う。  
【数５２】
2007509437.tif 000032
最小収束限界問題の目的関数Ｊ(ｕ(τ)，ｒ(τ)，Ｔ)の最小化により、Ｆ(ｙ(τ)，ｕ(τ)，ｒ(τ))を強制的に最小に収束する。この場合、最小の計算限界ＴにおいてＧ＝０であり、収束区間（convergence interval）は連続する反復変数の積分である。最適制御問題の解は、状態ベクトル変化（state-vector-change）の最適化の軌跡ｙ＊(τ)を生成するための最適制御ポリシーｕ＊(τ)およびｒ＊(τ)、τ∈[０，Ｔ]を提供する。最小収束限界問題の目的関数Ｊ(ｕ(τ)，ｒ(τ)，Ｔ)を定義する積分の最初の項によって、最適化の軌跡は連続型反復変数τに関してできるだけ効率的な軌跡ｙ＊に強制的に導かれる。以下の説明では、目的関数Ｆの偏導関数ベクトルＧとヘッシアン行列Ｈを頻繁に使用する。積分の第２項によって、ベクトルＧは強制的に０に向けられ、そのことによって目的関数Ｆは強制的に最小に導かれる。したがって、これらの項の組合せによって最適制御ポリシーｕ＊(τ)とｒ＊(τ)が生成され、このポリシーによって目的関数Ｆは最適な状態ベクトル変化の軌跡ｙ＊（τ）をたどって強制的に局所的最小値に導かれる。ここで、最適制御問題（１）〜（５）について説明する。第１の制約は正規化された状態ベクトル変化の最も急な降下を表しており、元の最適化問題が新しい最短時間制御問題に組み込まれる基本的な手段である。第２と第３の制約は、制御ベクトルの値ｕが常に非負であり、特定の極大制御値ｕｍａｘｉ未満であることと、ｕベクトルの値が最適化処理中に増大しないことを保証する。第４と第５の制約は、第２の内部的な制約と同様であり、束縛変数の係数ｒｉに関連する。第４と第５の制約は、束縛変数の係数ｒｉが最小値ｒｍｉｎｉを下回らないことと、束縛変数の係数ｒｉが最適化処理中に減少しないことを保証している。
【００４６】
最適制御理論の主な結果は、ポントリャーギンの最小原理（The Pontryagin Minimum Principle）として知られており、最小収束限界問題のフレームワークを提供し、状態ベクトル変化の最適化の軌跡ｙ＊（τ）を検出するための最適制御ポリシーｕ＊（τ）およびｒ＊（τ）、τ∈［０，Ｔ］を生成する。このアプローチには、ハミルトニアン（Hamiltonian）関数が含まれる。この関数は次のように定義される。  
【数５３】
2007509437.tif 000033
ハミルトニアンは、ｕ＊（τ）およびｒ＊（τ）におけるｕおよびｒの関数としての最小値をとる。換言すると、ハミルトニアンは最適制御ポリシーの最小値をとる。制御ポリシーｕ＊（τ）が最適な場合は、共役（costate）と呼ばれる絶対的に連続な関数Ｐ（τ）が存在する。このとき、状態関数ｙ(τ)およびと共役関数ｐ(τ)の導関数は、τ∈[０，Ｔ]となるすべてのτについて以下のように表され、
【数５４】
2007509437.tif 000034
以下の境界条件  
【数５５】
2007509437.tif 000035
を満たしている。制御問題の解は、次のように要約される。  
【数５６】
2007509437.tif 000036
【００４７】
このハミルトニアンは、機械系の総エネルギーを表現する機械系のハミルトニアンに類似している。関数ｙ＊（τ）は、連続型時間区間において機械系のコンポーネントのエネルギーが最小になる位置の指定に類似しており、共役関数ｐ＊（τ）は、連続型時間区間において機械系最小のコンポーネントのエネルギー・モーメントの指定に類似している。上で導かれた結果を使用すると、最適化の軌跡に沿った偏導関数∂ｙ＊／∂τの順方向（forward direction）の計算が可能である。偏導関数∂ｐ／∂τは、順方向の計算結果ｙ＊（τ）から、τ＝Ｔによって後退型で計算できる。もう一度伝統的な機械の問題に例えると、機械系のコンポーネントのモーメントは、本質的にコンポーネントの時間的な位置の変化によって決定する。したがって、モーメントまたは共役は、順方向で計算した位置から後退型で計算する必要がある。
【００４８】
ハミルトニアンの最小化は、それ自体が最適化問題であり、ハミルトニアンが最小化される場合の最適な制御ポリシーを生成する。最小化問題は、ニュートン近似法を使用して次のように解くこともできる。  
【数５７】
2007509437.tif 000037
したがって、最短時間問題および最適制御問題を定式化することにより、連続する反復変数τの０から限界Ｔまでの範囲にわたって、状態ベクトルｙおよび制御ベクトルｕとｒを順方向で計算し、共役ベクトルｐを後退型で計算する方法が得られる。さらに、この操作は１つの反復最適化法となり、ここで目的関数Ｆは、最適化制御ポリシーｕ＊（τ）およびｒ＊（τ）の下で最適な状態ベクトルの軌跡ｙ＊（τ）を経由して強制的に局所的最小値に導かれる。
【００４９】
図２０は、上記で導かれた状態ベクトル変化の最適化の軌跡ｙ＊および最適化制御ポリシーｕ＊とｒ＊に基づいた反復最適化法を図式的に説明している。反復最適化法は、２つのネストしたループ、すなわち反復変数Ｋで制御される外側ループと反復変数ｋで制御される内側ループからなると考えることができる。外側のネストしたループは、目的関数Ｆを最小化する状態ベクトルｙを検出するための、連続する繰り返しの範囲０からＴまでの状態ベクトル軌跡および制御ポリシーの偏微分の積分を表している。このような外側ループの繰り返しごとに、連続する反復変数τの短い区間に関して制御ポリシーベクトルｕとｒおよび状態ベクトルｙに関する複数の順方向の計算がステップτｋで実行される。ただし、ｋは０から最後のｋの値であり、換言すると、τ０、τ１、τ２、．．．、τｆｉｎａｌをとり、τｆｉｎａｌ≧Ｔになると、ｙ（τｋ）とｕ（τｋ）の計算結果を使用して、τ＝Ｔからτ＝０までの共役ベクトルｐ（τｋ）を後退型で計算する。内側ループを実行するごとに、ｙ（τｋ）の値の計算は前のｕ（τｋ）とｒ（τｋ）の計算に基づいて行われ、ｙ（τｋ）とｕ（τｋ）の両方の計算は直前の外側ループの繰り返しが終了したときに計算された共役ベクトルｐ（τｋ）を使用して行われる。最初の繰り返しでは、制御共役値は定数で近似される。
【００５０】
図２０において、外側ループは外側の円２００２で表され、内側ループの計算は、たとえばステップ２００４などの階段状の放射線と円周の断片で表される。繰り返しの最適化プロセスは、状態ベクトル、制御ベクトル、共役ベクトルの初期値ｙ０、ｕ０、ｒ０およびｐ０（それぞれ２００６）で開始される。外側ループ変数Ｋは０に設定される（２００８）。次に、複数の階段状の内側ループの繰り返しがｋ＝０で開始され、内側ループを次々に繰り返すたびに内側ループ変数ｋが１ずつ増分される。最初の内側ループの繰り返しで、制御ベクトルｕ（τｋ＋１）およびｒ（τｋ＋１）が計算される（２００４）。さらに、τ＝τｋ＋１、ｙ（τｋ＋１）における状態ベクトルの値が多くの連続型計算２０１０〜２０１６で計算される。最初の内側ループの繰り返しが終わるときに、局所的最小値に到達したかどうかが決定される。達した場合は反復最適化法が終了する。この様子は、図２０の破線２０１８で示されている。達していない場合は、次の内側ループの繰り返し２０２０〜２０２１が行われる。内側ループの繰り返しは、τｋ≧Ｔ（２０２２）となるまで続く。この時点で共役ベクトル値ｐ（τ）が、τ＝Ｔ  からτ＝０までの複数のステップ２０２４−２０３０で計算され、外側ループの繰り返し制御変数Ｋが増分されて完全な外側ループの繰り返しが終わり、次の外側ループの繰り返しを開始する点２０３２に戻る。
【００５１】
図１９と２０を併せて考察することによって、繰り返しによる最適化プロセスを図式的、概念的に見ることができる。最適化プロセスの内部ループの繰り返しは、そのそれぞれが図１９に示されている最適化の軌跡における一部分、たとえば一部分１９０６に対応すると考えることができる。外側ループ全体の１回の繰り返し（τ＝０からτ＝Ｔまで）によって、最適化の軌跡の大きな一部分が生成されると考えることができる（たとえば初期点１９０４から始まって点１９０８まで下方に伸びる最適化の軌跡の一部分）。したがって、複数の外側ループの繰り返しによって、複数の大きな最適化の軌跡の一部分を計算する。この一部分は、最終的に状態ベクトルｙを目的関数Ｆが少なくとも局所的最小値となる値に向かわせる。図２１は、反復最適化法全体を示す追加の図式的な説明を提供する。図２１に示されているように、初回の外側ループの繰り返しは、初期点２１０２から点２１０４までの最適化の軌跡の一部分を生成する。２回目の外側ループの繰り返しは、点２１０４から点２１０８までの最適化の軌跡の一部分２１０６を計算する。３回目の外側ループの繰り返しは、点２１０８から点２１１２に至る最適化の軌跡の一部分２１１０を生成する。限界の値Ｔは各外側ループの繰り返しごとに変化してもよいことに留意されたい。このように、図２１では、各最適化の軌跡一部分の終点はＴｋ（ただし、Ｋ＝１、２、および３）とマークされており、外側ループによる繰り返しの終点はそれぞれ次回の外側ループによる繰り返しの開始点となることが示されている。各内側ループのステップにおけるΔτのサイズも変化してよいことにも留意されたい。
【００５２】
図２２は、一般的な反復最適化のステップ（図２の２０３）を説明する制御流れ図である。まず、ｙ、ｕ、ｒおよびｐの初期値が選択され、外側ループの反復変数Ｋが初期値０に設定され、さらに初期限界Ｔ０がステップ２２０２で計算される。次に、内側ループを形成するステップ２２０４〜２２０８で、τ＝０からτｋずつ増分されてτ＝Ｔｋに至る間に、制御ベクトルｕとｒおよび状態ベクトルｙは順方向で伝達される。まず、ステップ２２０４で制御ベクトルｕとｒの次の値が計算される。さらに、ステップ２２０５で、τ＝τｋからτ＝τｋ＋Δτｋに状態ベクトルｙが伝達される。ステップ２２０６で最適化が収束したがどうかが判定される。収束した場合は、反復最適化が完了したので復帰する。収束していない場合は、ステップ２２０７でτｋ＋ΔτｋがＴＫより小さいかどうかが判定される。そうでない場合は、ステップ２２０８でτｋ＋ｌにτｋ＋Δτｋが代入され、内側ループの制御変数ｋが増分される。しかし、τｋがＴＫ以上の場合は、ステップ２２０９で共役の値がＴＫから０まで逆順に計算され、ステップ２２１０で、反復変数および状態ベクトル、制御ベクトルの計算された値を含むデータベースが更新される。
【００５３】
凸問題（convex problems）では臨界点は１つしかないので、図１９〜２１を参照しながら以上で説明した反復最適化法は、数千もの決定変数を持つ極めて大きな超次元問題領域に対しても、通常は妥当な時間内に唯一の軌道で収束する。しかし、適用例における最適化問題は通常凸ではなく、問題領域は多くの臨界点を持つという特徴がある。図２３は、このような問題領域を示している。図２３において、問題領域２３０３の境界面は、複数の中央が突出した突起部（たとえば突起部２３０４）を持つ椀状の形である。これらの突起部は、３つの頂点（２３０４、２３０６、２３０８）がある小さな山岳地帯と、深い局所的な谷、すなわちくぼみ２３１０と考えることができる。問題領域の大域的な最適化の価値は、外側のくぼみの真中にある山岳地帯の下方深くに存在する。
【００５４】
図２４は、図２３と同じ問題領域を示しているが、臨界点が明示的に示されている。山岳地帯２３０４、２３０６および２３０８の頂上は問題領域内の極大を表しており、小さな谷の底２３１０は局所的最小値を表している。山岳地帯内の鞍点（たとえば鞍点２３１２および２３１４）も臨界点である。図２５は、臨界点（たとえば図２４における臨界点２３０４）の拡大表示である。図２５に示すように、臨界点２５０６に近いいくつかの点２５０２〜２５０４で計算されたｘ方向およびｙ方向の偏微分はすべて絶対値が小さい。したがって、臨界点２５０６に近い点２２０２〜２３０４で計算された勾配も０に近い。要約すると、τに関する目的関数の導関数
【数５８】
2007509437.tif 000038
によって、臨界点に近い最適化の軌跡を決定する手段が得られる。したがって、点が臨界点に近づくにつれて勾配が０になるかまたは０に近づくと、τに関する状態ベクトルの軌跡の変化（ｒａｔｅ）も０となる。さらに、Ｇは臨界点の近くでは敏感であり、したがって小さな局所的な変化がＧの比較的大きな変化の原因となることがあり、Ｇの符号が変化する場合もある。しかし、Ｇが小さい場合は、状態ベクトルの軌跡を計算するために使用するニュートン法や他の同様の方法は、一般的に最適化の軌跡を発見するためには有効ではない。換言すると、勾配フィールドは一般的に準最適軌跡または最適な軌跡の方向を提供する。したがって、勾配フィールドが０の場合は、可能なすべての軌跡が同等に適切であるように思われる。これは、磁北極または磁南極でコンパスを使おうとする場合と同様の状況である。後で説明するように、臨界点の近くでは別の方法を使用して状態ベクトルの最適化の軌跡を計算する必要がある。
【００５５】
臨界点については、別の考慮すべき問題がある。この考慮すべき問題には、問題領域におけるいくつかの可能性のある局所的最小値の識別と評価を試みるために、大域的な最適化法で臨界点の近くで使用する必要がある検索に似た戦略が含まれる。図２６は、第１の検索技術を示している。図２６では、水平線Ｔｎへ向かう本最適化法の外側ループの繰り返しは状態ベクトルを点２６０２に伝達する。この点は、山岳地帯の頂上である臨界点２３０６の近くである。この点における勾配が０に近いという事実によって、最適化法に対して臨界点の近くに到達したことを伝える。ただし、次のステップは到着した臨界点のタイプによって決まる。たとえば点２３０６のように極大の場合は、本最適化法では何らかの方法を使用して最適化の軌跡を大きく変更し、極大２３０６を離れて局所的最小値に向かう降下２３０８を続ける必要がある。しかし、たとえば局所（local）２３１０のような局所的最小値の近傍では、本最適化法は局所的最小値２３１２から別の点２３１４に上昇し、この点から局所的最小値に向かって降下２３１６を続けるような技術を使用する必要がある。
【００５６】
図２７は、問題領域を検索するための第２の戦略を示している。図２７では、検索方法は極大２３０６に近い点２６０２に再度到着している。しかし、極大を検出すると、本最適化法は異なる方向２７０２〜２７０５に向かって降下を続ける相対的に独立した最適化の軌道をいくつか追加して生成する。相対的に独立した軌道の１つである２７０４は、鞍点２３１４に近い点２７０８に到達する。この点で、この軌道２７０４は追加の軌道を分岐または生成し、結果として得られる相対的に独立した各軌道は鞍点から異なる方向２７１０と２７１２に向かって降下する。最初に生成された別の相対的に独立した軌道２７０２は、局所的最小値２３１０に到達する。この時点で、新しい最適化の軌道を生成せず、相対的に独立した軌道２７０２は終了し、親軌道２７０１に局所的最小値を報告する。親軌道は子によって識別された最小の局所的最小値を選択し、それを局所的最小値としてその親に返す。結果として問題領域の再帰的な検索が行われ、最適化問題の準最適解または最適解が識別される。別の実施形態では、親軌道は報告された局所的最小値の情報を使用して問題領域を切り詰め、局所的最小値に基づいて問題領域を修正でき、問題領域が切り詰められたことを相対的に独立した他のすべての軌道に伝えることができる。
【００５７】
本最適化法は、図２６に示す唯一の軌道の技術でも、図２７に示す複数の軌道の技術でも、勾配が０に近い値に減少したときに到達した臨界点の種類を決定する手段が必要である。臨界点の種類を解明するための臨界点指標（critical-point beacon）は前述のヘッシアン行列である。勾配が０に近づいたときに、ヘッシアン行列のすべての要素が負ならば、極大に到達している。勾配が０で、ヘッシアン行列のすべての要素が正ならば、局所的最小値に到達している。それ以外の場合に、勾配が０でヘッシアン行列の要素の符号が混在していれば、変曲点または鞍点に到達している。正則（non-singula）ヘッシアン行列の臨界点指標属性は形式的に次のように記述される。
【００５８】
点ｙについて、
【数５９】
2007509437.tif 000039
のとき、
Ｈの全ての固有値が負ならば、ｙは局所的または大域的な最大値の近傍にあり、
Ｈの全ての固有値が正ならば、ｙは局所的または大域的な最小値の近傍にあり、
その他の場合、ｙは鞍点（saddle point）にある。
【００５９】
図２６と２７に関連して前述したように、唯一の最適化の軌道は、階層的に体系化されたモジュールの集合で構成される。このようなモジュールを離散型の計算環境で同時に実行し、前述の最適化法による状態ベクトルの伝達に関連する大量の計算タスクを分散してもよい。図２８は、１つの最適化の軌道を構成するモジュールとモジュール化された編成を示している。本明細書では、最上位レベルのモジュールを「スーパー（super）」２８０２と呼ぶ。スーパーには、中間結果の比較的大きなデータベース２８０４が含まれる。スーパーは、それぞれの連続する限界ＴＫの計算、図２７を参照しながらすでに説明したように、初期点の選択（複数の軌道による最適化法におけるルート軌道の場合）、および下位モジュールの起動とアクティビティ調整を担当する。第１レベルの下位モジュールには、ｕ＿ａｎｄ＿ｒモジュール２８０６、ｙモジュール２８０８およびｐモジュール２８１０が含まれる。ｕ＿ａｎｄ＿ｒモジュール２８０６は、図２２のステップ２２０４に示すように、制御ベクトルｕとｒの順方向の計算を担当する。
【００６０】
ｙモジュール２８０８は、図２２のステップ２２０５に示すように、連続する区間Δτを通じた状態ベクトルの伝達を担当する。ｐモジュール２８１０は、図２２のステップ２２０９に示すように、反復最適化法の外側ループの繰り返しが終了するたびに、Ｔｋから０までの後退型による共役ベクトル値の計算を担当する。
【００６１】
ｕ＿ａｎｄ＿ｒモジュール２８０６は、制御ベクトルの順方向の計算を、ｍ個のエージェント２８１２〜２８１５に分散している。ｙモジュール２８０８とｐモジュール２８１０は各モジュールのタスクをそれぞれｎ個のエージェント２８１６〜２８１８および２８２０〜２８２３に分散している。ｙモジュール２８０８がｙの伝達をｎ個のｙモジュール・エージェント２８１６〜２８１８に分散することは、最適化計算全体を分散し、パラレルな計算の機会を提供するために重要なアプローチである。ｙモジュール２８０８およびそのエージェント２８１６〜２８１８で行われる状態ベクトル伝達の計算は重要であり、最も計算コストがかかり、状態ベクトル変化による最適化の軌跡の計算に必要なタスクである。本発明の別の実施形態では、ｐモジュールは効率的な補間法を使用して共役ベクトルの値を推定するので、問題をｎ個のｐモジュール・エージェントに分散する必要はない。したがって、ｐモジュールは分散せずに共役を計算する。さらに別の実施形態では、ｐモジュールは共役の補間を少数のｐモジュール・エージェントに分散できる。
【００６２】
スーパー２８０２で保守するデータベース２８０４は大きいため、メモリだけでなく、軌道のモジュールを分散する１つまたは複数のコンピュータ・システムに直接接続またはリモート接続されたマス・ストレージ・デバイスに格納してもよい。図２９は、本発明の様々な実施形態においてスーパーが保守するデータベースの内容の一部を示している。一般に、スーパーは、０からＴまでの水平線の各部分区間τｋで、モジュールによって計算されるさまざまな量の値を保守する必要がある。したがって、図２９では、たとえばベクトルＧのデータ構造２９０４の部分列（sub-column）２９０２など、データ構造の第１列の各部分列は、特定のτの値に対応するエンティティの計算値を表している。図２９において、部分列２９０２はτ＝０におけるベクトルＧの値を表している。スーパーは、状態ベクトルの値２９０６、制御ベクトルｕとｒの値２９０８〜２９０９、ベクトルｗの値２９１０、スカラーρの値２９１１、モース（Morse）フラグ・ベクトルの値２９１２、共役ベクトルｐの値２９１４、目的関数Ｆのスカラー値２９１５、ベクトルＧの値２９０４、ハミルトニアンの値２９１６、現在のＴと反復制御変数および他のそのような局所的な制御変数の値２９１８、およびヘッシアン行列の逆行列の値２９２０〜２９２７を保存する。
【００６３】
図３０〜３４は、図２８の描画法を使用して、本最適化法の外側ループが繰り返されるときに、階層的に構造化された最適化の軌道のモジュール間で情報が流れる様子を示している。ルート軌道の場合は、図３０に示すように、最適化の軌跡の初期点が決定する（３００２）。さらに、他のデータ構造の中で、状態ベクトル、共役ベクトル、制御ベクトル、および制御変数の初期値が決定され、さまざまなデータ構造と要素３００４が初期化される。ルート以外の軌道は初期点および、状態ベクトルや制御ベクトルなどのさまざまな計算エンティティの現在の値を親軌道から継承する。図３１に示すように、内側ループの最初のステップ（図２２の２２０４）で、ｕ＿ａｎｄ＿ｒモジュール２８０６はスーパー  ３１０２から状態ベクトル、共役ベクトル、ベクトルｗおよびヘッシアンの逆行列の値を受け取り、これらの値および／またはこれらの値の一部をｕ＿ａｎｄ＿ｒモジュール・エージェント３１０４〜３１０７に分散し、中間の計算結果をｕ＿ａｎｄ＿ｒモジュール・エージェント３１０８、３１０９、３１１０、および３１１１から受け取り、新しい制御ベクトルの値を計算し、さらに新しい制御ベクトルの値３１１２をスーパーに返す。図３２は、内側ループの次のステップ（図２２の２２０５）におけるデータの流れを説明している。このステップでは、ｙモジュール２８０８はスーパー  ３２０２から新しく計算された制御ベクトルｕとｒを受け取り、これらの値にヘッシアンや他のベクトルなどの値を含むスーパーから取得したその他の値および／またはｙモジュールが計算したその他の値を追加してｙエージェントに配信する（３２０４〜３２０６）。次に、中間の計算値をｙエージェントから受け取り（３２０８、３２０９、３２１０）、特に、状態ベクトル、ｗベクトル、ｗベクトルの状態ベクトルに関する微分、ヘッシアンの逆行列を計算してスーパーに返す（３２１０）。ｙモジュールは、前述のように、内側ループの繰り返しのたびに、多くの個別の伝達ステップを実行する。図３３は、図３１と３２に示すデータ交換ステップが内部ループの反復回数だけ続けて繰り返される様子を示している。図３４は、図２２のステップ２２０９で、ｐモジュール２８１０で実行する共役ベクトルの計算に関連するデータ交換を示している。ｐモジュールは、スーパー  ３４０２から状態ベクトル、制御ベクトルおよび他の計算値を受け取り、これらの値またはこれらの値の一部をｐエージェントに配信し（３４０４〜３４０７）、中間結果をｐエージェントから受け取り（３４０８〜３４１１）、τの値０からＴｋまでの範囲の共役ベクトルを計算してスーパーに返す（値３４１２）。前述のように、現在の好ましい実施形態では、ｐモジュールは共役ベクトルを積分せずに補間するので、図３４に示すような方法で問題を分散する必要はない。
【００６４】
以下の説明では、表記上の都合により、ベクトルｕとｒを併せてベクトルｕとする。つまり、後続の説明では、束縛変数の係数とバリア関数の係数が統合され、唯一の制御ベクトルｕを構成する。本発明の実施形態で採用する反復最適化の戦略には、
【数６０】
2007509437.tif 000040
の積分が含まれることを想起されたい。ここで、
ｗ(ｙ(τ)，ｕ(τ))＝(Ｈ(ｙ(τ)，ｕ(τ)))−１Ｇ(ｙ(τ)，ｕ(τ))
ｐ(ｙ(τ)，ｕ(τ))＝(Ｇ(ｙ(τ)，ｕ(τ)))Ｔｗ(ｙ(τ)，ｕ(τ))
である。一実施形態において、ｄｙ(τ)／ｄτの積分は、順方向オイラー積分法（forward Euler integration method）によって可変ステップサイズΔτｋで、以下のように実行される。  
【数６１】
2007509437.tif 000041
したがって、上の式で明らかなように、ｙモジュールは順方向積分のステップごとにヘッシアン行列の逆行列を求め、ヘッシアンの逆行列とベクトルＧの積としてｗベクトルを計算する必要がある。図３５は、この計算を図式的に示している。ヘッシアン行列の逆行列は、最適化に必要なタスクの中で計算量が群を抜いて最大である。したがって、ｙモジュールはヘッシアン行列の逆行列の計算をｙモジュール・エージェントに分散する。ただし、この分散処理はいくつかの仮定に基づいている。
【００６５】
以下の式を想起されたい。  
Ｈ(ｙ(τｋ)，ｕ(τｋ))ｗ(ｙ(τｋ)，ｕ(τｋ))＝Ｇ(ｙ(τｋ)，ｕ(τｋ))
図３６は、ヘッシアン行列とｗベクトルの乗算が分散方式でｇベクトルを生成する様子を示している。図３６に示すように、Ｇベクトルとｗベクトル、およびヘッシアン行列は区分ｉ＝０から５まで縦方向に分割できる。Ｇベクトル３６０２の最初の区分をＧ０とし、ヘッシアン行列３６０４の最初の縦区分をＨ０とし、さらにｗベクトル３６０６の最初の区分をｗ０とする。ヘッシアン行列は、さらに横方向に横区分ｊ＝０から５に分割できる。図３６において、ヘッシアン行列３６０８の最初の横区分はヘッシアン行列内の最初の縦列として表示されている。
【００６６】
ヘッシアン行列を縦と横の両方に区分化すると、ヘッシアン行列がブロック状の部分行列に分割される。たとえば、図３６で、ヘッシアン行列が左上の部分行列Ｈ０，０（３６１０）から右下の部分行列Ｈ５，５（３６１２）までの３６個の部分行列に区分化されることを示している。対角線上のブロック部分行列Ｈ０，０，Ｈ１，１，Ｈ２，２，Ｈ３，３，Ｈ４，４、およびＨ５，５が他の対角線上にない部分行列の要素よりもはるかに絶対値の大きな要素を持つ対角優位（diagonally dominant）ヘッシアン行列の場合は、ヘッシアン行列とｗベクトルとの行列乗算をｗベクトルの区分の行列とヘッシアン行列の対応する縦区分内にある対角線ブロック部分行列との乗算で近似してもよい。たとえば、行列全体の乗算を行う場合は、Ｇベクトルの区分Ｇ０はｗベクトルの区分とヘッシアン行列の最初の縦区分の対応するブロックとの行列乗算による積の和に等しい。これは、次のように表される。  
【数６２】
2007509437.tif 000042
【００６７】
ヘッシアン行列が対角優位である場合は、対角上にない部分行列は要素の絶対値が相対的に小さいため、ｗベクトル区分の計算における重要度は低い。この状況は、ヘッシアン行列が各区分と他すべての区分との間にある依存関係をできるだけ多く削除する方法で区分化された場合に発生する。したがって、特定のｙモジュール・エージェントは、対応するヘッシアン行列の対角ブロック部分行列の逆行列３８０８を計算し、大きなカッコでくくられた項（図３８の３８１２）をブロック対角線ヘッシアン部分行列の逆行列３８０８にかけることによって、ｗベクトル区分３８１０の現在の値を近似して計算できる。ただし、残りのｗベクトル区分３８１４〜３８１７については前に計算された値を使用する。
【００６８】
図３９は第２のｙモジュール・エージェントによるｗベクトル区分ｗ１の計算を説明する。第２のｙモジュール・エージェント３９０２は、現在のヘッシアン行列の対角ブロック部分行列Ｈ１，１と現在のＧベクトル区分Ｇ１および最新に計算されたｗベクトル区分ｗ０、ｗ２、ｗ３、ｗ４およびｗ５をｙモジュール２８０８から取得し、現在のヘッシアンブロック対角部分行列Ｈ１，１の逆行列を計算してから、計算されたヘッシアン行列の逆行列を使用してｗベクトル区分ｗ１の現在の近似値を計算する。次の計算ステップで、ｙモジュール・エージェント３９０２は、ｗベクトル区分の現在の近似値を他のｙモジュール・エージェントに配信できるようにｙモジュール２８０８に提供でき、次にτが増分されたときに次のｗベクトル区分ｗ１の値を計算できるように、他のｙモジュール・エージェントが最新に計算した他のｗベクトルの区分ｗ０、ｗ２、ｗ３、ｗ４およびｗ５の値を受け取ることができる。この方法を使用すると、ｙモジュール・エージェントは各ステップで同期を取りながらｗベクトル区分の近似値を計算できる。
【００６９】
次に、スーパーモジュールの制御流れ図について説明する。図４０はスーパーモジュールの最も高い概念レベルを示す制御流れ図である。ステップ４００２で、スーパーモジュールは、自らが親軌道から生成された複製または親軌道から分岐した子軌道であるかどうかを判断する。そうでない場合、したがって現在の軌道がルート軌道である場合は、スーパーモジュールはステップ４００４で最適化の初期点を決定する。ステップ４００６で、状態ベクトル、共役ベクトル、および他の値のセットを含むさまざまな計算エンティティを初期化し、さらにブール変数「critical」および「precritical」のブール値ＦＡＬＳＥへの初期化を含む多くの変数の初期化も行う。ルート軌道による初期化に続いて、あるいは複製された軌道の場合にステップ４００２からの制御の流れに従って進んだ結果として、スーパーモジュールは次にステップ４００８で限界Ｔを選択する。さらに、スーパーモジュールは前述の最適化法における連続する内側ループの繰返しを実行する関数「ｕおよびｙをＴまで積分（integrate_u_and_y_to_T）」を呼び出す（ステップ４０１０）。次にステップ４０１２に進み、スーパーモジュールはステップ４０１０のルーチン「ｕおよびｙをＴまで積分」で実行した内側ループの繰り返しに基づいて決定を行うために、ルーチン「評価（evaluate）」を呼び出す。ルーチン「評価」は、前述の最適化法の外側ループの繰り返しをさらにもう１度実行するかどうかを決定するために、ステップ４０１４でスーパーモジュールが評価するブール値「継続（continue）」を返す。実行する場合は、スーパーモジュールは、計算した共役ベクトルの値を返すためにルーチン「ｐを積分（integrate_p）」４０１６を呼び出してから、ステップ４００８に戻って新しい外側ループの繰り返しを開始する。スーパーモジュールが外側ループの繰り返しを継続しないと決定された場合は、スーパーモジュールはステップ４０１８でルーチン「報告（report）」を呼び出してから復帰する。
【００７０】
次に、ステップ４００４、４００８、４０１０、４０１２および４０１８を示す制御流れ図について説明する。ステップ４０１６で呼び出されるルーチン「ｐを積分」と、ルーチン「ｕおよびｙをＴまで積分」から呼び出されるルーチン「ｕを積分（integrate_u）」および「ｙを積分（integrate_y）」は、明確かつ簡潔にするために、制御流れ図を使用せずに数理的に記述する。
【００７１】
図４１は、図４０のステップ４００４で呼び出されるルーチン「初期点（initial_point）」の制御流れ図である。初期点ルーチンは、状態ベクトルｙ０の初期値を決定するために、ステップ４１０２でユーザーから、あるいはオプティマイザの数理モデル作成からまたは前処理の部分から提供されるすべてのヒントまたは経験則を使用して値を推測する。次に、初期点ルーチンは、ステップ４１０４〜４１０７で構成されるｆｏｒループの中で、各制約ｇｉを評価し、最初に推測した状態ベクトルの点ｙ０で評価された制約が０より大きいかどうかを判定する。大きい場合は、ステップ４１０７で区分Ｉに制約が指定される。大きくない場合は、区分Ｊにステップ４１０６で制約が指定される。さらにステップ４１１０で、初期点ルーチンは目的関数を次のように構成する。  
【数６３】
2007509437.tif 000043
これは、ステップ４１１２でオプティマイザを再帰的に呼び出し（図３０）、修正された初期状態ベクトルの値を生成することによって最小化される。目的関数Ｆ（上部にチルダ）は、初めは状態ベクトルの初期値を包み込むために目的関数Ｆの問題領域境界を拡張すると考えることができ、次にＦ（チルダ）を最適化することによって問題領域境界を縮小し、初期点を強制的に目的関数Ｆの問題領域の内部に含める。ステップ４１１４で、初期点ルーチンは変更した状態ベクトル値が実現可能かどうかを判断する。実現可能な状態ベクトル値は、局所的最小値を検出する最適化法で局所的最小値に十分近づく可能性の高い問題領域内の点である。変更した状態ベクトル値が実現可能な場合は、最適化の初期点として返される。それ以外の場合は、制御フローはステップ４１０２に戻り、状態ベクトルの初期値を決定する他のパス（pass）を検出する。後続の繰り返しでは、ステップ４１０２で前の繰り返しから収集した情報を使用して、状態ベクトルの初期値のよりよい推定値を得る。
【００７２】
図４２は、ルーチン「次のＴを選択（choose_next_T）」の制御流れ図である。ステップ４２０１では、ルーチン「次のＴを選択」は最小の限界定数に局所変数「Ｔ」を割り当てる。次にステップ４２０３〜４２０７で構成されるｆｏｒループで、ルーチン「次のＴを選択」は、ｙモジュール・エージェントに関連したヘッシアン対角ブロック部分行列の最小絶対固有値で定数パラメータを割って得られたＴの値をｙモジュール・エージェントごとに計算する。ｙモジュール・エージェントについて計算したＴの値がステップ４２０５で検出した変数「Ｔ」の内容より大きい場合は、変数「Ｔ」にｙモジュール・エージェントについて計算されたＴの値が割り当てられる。次にステップ４２０８で、ルーチン「次のＴを選択」は局所変数「ｎ」を増分の最大数に初期化し、局所変数「ａ」を増分係数に初期化し、さらにループ制御変数「ｊ」を０に初期化する。ステップ４２１０〜４２１２で、ルーチン「次のＴを選択」は、変数「Ｔ」の現在の内容に対して評価されたＴに関する共役ベクトルの微分が０に近づくまで（ステップ４２１０で判定する）、ステップ４２１２で変数「Ｔ」の値を増分する。このループは、繰り返しの回数が局所変数「ｎ」の内容に等しくなった場合にも終了する。ｄｐ／ｄＴ＝０の場合は、共役ベクトルの計算を０〜Ｔの範囲で継続して繰り返すことによって、より容易に計算を実行できる。本発明の代わりの実施形態では、複雑な積分を使用せずに、多項式による近似を使用して共役ベクトルを計算できる。
【００７３】
図４３は、図４０のステップ４０１２で呼び出されたルーチン「評価」の制御流れ図である。ステップ４３０２で、「評価」のルーチンはブールフラグ「at_minimum」および「continue」にブール値ＦＡＬＳＥを設定する。次にステップ４３０４で、ルーチン評価は最適化に対して割り当てられたコンピュータ・リソースが使い果たされたかどうかを判断する。これには、指定された最大時間より長い時間継続する計算、処理中のリソースの使用、メモリ、または使用する指定の最大リソース量よりも多いコンピュータ・リソースの組合せ、エラー状況の発生、または他の種類のあらかじめ指定された終了基準を含めてもよい。コンピュータ・リソースが使い果たされている場合はルーチン評価が復帰する。それ以外の場合は、ステップ４３０６でルーチン評価は最適化の軌道が現在臨界的（critical）であるかどうかを判断する。すべてのｙモジュール・エージェントが現在モースモードであり、あらかじめ指定された回数の繰り返しの間そのモードであった場合、最適化の軌道は臨界的である。モースモードについては、以下でより詳細に説明する。現在の軌道が臨界的でない場合に、ルーチン「評価」はステップ４３０８でブールフラグ「継続」をＴＲＵＥに設定して返す。臨界的である場合は、ステップ４３１０で、ルーチン「評価」は現在臨界的であると判断されている軌道が局所的最小値に到達しているかどうかを判断する。到達している場合は、ステップ４３１２で、ブールフラグ「ａｔ＿ｍｉｎｉｍｕｍ」にＴＲＵＥが設定され、ルーチン「評価」は復帰する。最適化の軌道が局所的最小値ではないが臨界的である場合は、ルーチン「評価」はステップ４３１４で生成する複製軌道の数と、複製軌道の最初の方向または開始点なども決定する。複製軌道の生成は、図２７に示されている。ステップ４３１６〜４３１９のｆｏｒループでは、指定された数の複製軌道が生成され、リンクされた各複製軌道には、現在の最適化の軌道の子が存在する。ステップ４３２０で、ルーチン「評価」は現在の最適化の軌跡を摂動させ（perturb）、この臨界点から最適化の軌跡を離して局所的最小値に向けて降下を続けるようにしてもよい。ステップ４３２２で、ルーチン「評価」はブールフラグ「継続」にブール値ＴＲＵＥを設定してから復帰する。
【００７４】
図４４は、ルーチン「報告」の制御流れ図である。このルーチンは、図４０のステップ４０１８で呼び出される。ステップ４４０２で、ルーチン「報告」は現在の最適化軌道に親があるかどうかを判断する。親がない場合は、現在の最適化の軌道がルート最適化軌道であれば、ステップ４４０４でルーチン「報告」は現在のルート最適化軌道にリンクしているまだ動作中の子軌道があるかどうかを判断する。子がある場合は、現在のルート最適化軌道はステップ４４０６で子最適化軌道が実行を終了するまで待機する。次に、ステップ４４０８で、ルーチン「報告」は現在の状態ベクトル値で計算した目的関数の現在の値を局所変数「ｍｉｎ」に設定し、局所変数「location」に現在の状態ベクトルの値を設定し、局所変数「at_min」に現在の最適化軌道に関してルーチン「評価」から返されたブール値フラグ「at_minimum」を設定する。ステップ４４１０で、ルーチン「報告」は、子軌道が終了し、現在の最適化の軌道に報告したかどうかを判断する。報告した場合は、ステップ４４１２〜４４１５のｆｏｒループで、現在の最適化の軌道は任意の子から報告され、終了した各子が生成したレポートを評価して、現在実行中の最適化の軌道により達成された目的関数に対する最小の値を決定する。次にステップ４４１６で、ルーチン「報告」は実行中の子がまだあるかどうかを判断する。ある場合は、ステップ４４１８〜４４２０のｆｏｒループで、ルーチン「報告」は各子を現在実行中の軌道の親にリンクする。ステップ４４１６のテストは、現在実行中の最適化の軌道がルート軌道である場合に失敗することに留意されたい。ステップ４４２２で、現在実行中の最適化の軌道がルート軌道であると決定された場合は、ステップ４４２４でルーチン「報告」は最小の目的関数値、状態ベクトル、およびこの最適化が現在最小であるかどうかの指示を返す。ルート以外の場合は、ルーチン「報告」はステップ４４２６で前述の同じ値を親軌道に報告する。
【００７５】
図４５は、図４０のステップ４０１０で呼び出されるルーチン「ｕおよびｙをＴまで積分」の制御流れ図である。ステップ４５０２〜４５０７のｆｏｒループでは、最適化法による内部ループの１回の実行を完了するために、ルーチン「ｕおよびｙをＴまで積分」はステップ４５０３と４５０４でルーチン「ｕを積分」と「ｙを積分」を連続して呼び出す。連続型繰り返しは、ステップ４５０５で最適化の軌道が臨界的になったか、値τｋがＴより大きいか等しくなったと判定されるまで、またはステップ４５０７で、コンピュータ・リソースの不足などその他の停止基準が検出されるまで行われる。
【００７６】
図４６は、ルーチン「ｙを積分」の制御流れ図である。このルーチンはｙモジュール（図２８の２８０８）の動作を説明する。ステップ４６０２で、ルーチン「ｙを積分」はルーチン「ｙを積分」が必要とするヘッシアン行列、ベクトルＧ、およびその他のデータのシンボリック表現をスーパーモジュールから取得する。次にステップ４６０４で、ルーチン「ｙを積分」はヘッシアン行列とベクトルＧの数値を計算する。ステップ４６０６では、ルーチン「ｙを積分」は各ｙモジュール・エージェントへの計算済みのヘッシアン行列とＧベクトルの区分を配信する。ｙモジュール・エージェントは、関連するｗベクトル区分の値と配信された行の値を計算し、ステップ４６０８でルーチン「ｙを積分」がこれらを収集する。ステップ４６１０で、ルーチン「ｙを積分」はエージェントが現在モースモードであるかどうかを判断する。
【００７７】
用語「モースモード（Morse mode）」は、現在エージェントがベクトルＧのエージェントの区分に関連する臨界点の近くにあるので、現在ｙエージェントの方法に関連する当初の目的関数のモース二次近似に基づいてそのｗベクトルの区分を計算する必要があるエージェントを表す。いずれかのエージェントがモースモードである場合、ステップ４６１２ではルーチン「ｙを積分」はモースモードのエージェントｉに対して、ρｉおよびｗｉの再計算を指示し、再計算したρｉとｗｉをモースモードの各エージェントｉから収集する。次にステップ４６１４で、ルーチン「ｙを積分」はすべてのエージェントが現在モースモードであるかどうかを判断する。すべてのエージェントがモースモードの場合は、ステップ４６１６でルーチン「ｙを積分」は局所変数「precritical」の現在の値がブール値ＴＲＵＥかどうかを判断する。ＴＲＵＥの場合は、エージェントは以前に実行された繰り返しですでにモースモードであるので、ステップ４６１８で局所変数「critical_count」を増分する。ステップ４６１９で、ルーチン「ｙを積分」は局所変数「critical_count」の内容をしきい値定数と比較し、局所変数「critical_count」の内容が臨界的な定数と等しいか超えている場合は、ルーチン「ｙを積分」はステップ４６２０でブールフラグ「critical」にブール値ＴＲＵＥを設定することによって、現在の最適化の軌道が臨界的であることを宣言するのに十分な回数の繰り返しにわたって臨界的な状態であったと結論付ける。しかし、変数「critical_count」の内容が現在の最適化の軌道が臨界的であると宣言するための正当な理由にはならない場合は、ステップ４６２４に制御を渡す。ステップ４６１６に戻ると、ルーチン「ｙを積分」が変数「precritical」の値がブール値ＴＲＵＥではないと判断した場合、ステップ４６２２では変数「precritical」は真に設定され、変数「critical_count」は０に初期化されてから、制御はステップ４６２４に移る。ステップ４６２４では、ルーチン「ｙを積分」は広域的なρとｗベクトル値を計算する。それらの値からｄｙ／ｄτの値を計算する。
【００７８】
次に、ステップ４６２５で、ルーチン「ｙを積分」は現在のオイラー順方向積分（Euler forward integration）のステップに関してΔτを決定する。Δτの計算が成功すると（ステップ４６２６で決定する）、ルーチン「ｙを積分」は前進積分ステップを実行し（ステップ４６２８）、必要に応じて束縛変数の係数を更新して（ステップ４６３０）復帰する。Δτの計算が失敗すると（ステップ４６２６で決定する）、ルーチン「ｙを積分」はエラーで復帰する。
【００７９】
オイラー順方向積分の技術について前述したが、形式的には次のように表される。ρの広域的な値は、現在の内部ループ・ステップでｙモジュール・エージェントによって計算された区分化されたρの値から、次のように組み立てられる。  
【数６４】
2007509437.tif 000044
次に、ｙモジュール・エージェントによってｗベクトルの区分の値が収集され、組み立てられてｗベクトルの値が生成される。さらに、後述する方法で連続型反復の増分Δτが計算される。最後に、状態ベクトルは現在の連続型反復の値τｋから次の連続型反復の値τｋ＋Δτに次のようにして伝搬される。  
【数６５】
2007509437.tif 000045
最終的に、連続型反復の値τｋは次のように更新される。  
τｋ＝τｋ＋Δτｋ
ただし、ｙモジュール・エージェントの局所的なρの絶対値であるρｉが特定の定数パラメータを下回る場合は、ｙモジュール・エージェントは、ｙモジュールのベクトルＧ区分に関して臨界的とみなされることに留意されたい。ｗベクトルの区分ｗｉとｙモジュール・エージェントの局所的なρであるρｉの計算については以下で説明する。
【００８０】
図４７は、図４６のステップ４６２５で呼び出すルーチン「Δτを計算（compute Δτ）」を示す制御流れ図である。ステップ４７０２で、Δτの見積もりまたは推定値が決定する。この見積もりまたは推定値は、定数パラメータから求めることができる。または、より複雑な実施形態では、以前に計算したΔτの値と状態ベクトルの値、制御ベクトルの値、および他の同様の情報の最新の履歴から、経験的な方法（empirical methods）または経験則（heuristics）によって導くこともできる。説明されている実施形態において、オイラー順方向積分のステップで制御ベクトルがほぼ一定の値をとる場合は、次のようにＦがステップごとに１ずつ減少する。  
【数６６】
2007509437.tif 000046
ｙモジュール・エージェントのいくつかがモースモードの場合は、ｄＦ／ｄτが値「−１」を外れてもよい。しかし、ｄＦ／ｄτの値が０を上回る場合は、ルーチン「Δτを計算」のステップ４７０４で検出するように、エラー状態が発生する。これは、ｄＦ(τｋ)／ｄτが降下の方向ではないので、オイラー順方向積分では実現可能なｙ＊（τ）の軌跡を生成できないためである。それ以外の場合は、Δτの推定値Δτ'がオイラー順方向積分のステップ式、  
【数６７】
2007509437.tif 000047
に挿入され、状態ベクトルがｙ（τｋ＋Δτ'）に伝搬される。この計算された次の状態ベクトルの値についてその実現可能性がチェックされる（ステップ４７０６）。状態ベクトルの値が実現可能でない場合は、Δτの推定値Δτ'が減少し（ステップ４７０８）、新しい前進状態ベクトルの値が再計算され、その実現可能性がチェックされる（ステップ４７０６）。実現可能な前進状態ベクトルと対応するΔτの推定値Δτ'が検出された場合は、ステップ４７１０でArmigoの条件、
【数６８】
2007509437.tif 000048
がΔτの推定値Δτ'に関連して評価される。Ａｒｍｉｇｏの条件が満足された場合は、Δτの推定値Δτ'が新しく計算されたΔτとして返される。それ以外の場合は、Δτが近似される。複数の異なる近似方法を利用できる。Δτの以前の実現可能な推定値ΔτＡがある場合は、次のようにτｋを中心とする目的関数Ｆの２次近似を作成できる。  
【数６９】
2007509437.tif 000049
【００８１】
次に、Δτの推定値が次のようにして得られる。  
【数７０】
2007509437.tif 000050
Ｆを最小化する降下の方向に沿ったΔτは、Ｆ（チルダ）を最小化するΔτを検出することによって近似できる。これは次の場合に発生する。  
【数７１】
2007509437.tif 000051
したがって、次の式が得られる。
Δτ’＝―ｂ／２ａ
現在のτｋで、以前に計算されたΔτの２つの実現可能な推定値が既知の場合は、３次近似法を利用してもよい。ただし、ステップ４７１２でΔτが近似され、新しい近似値Δτを使用して、状態ベクトルの前進した値が再計算されると、ステップ４７０６でその実現可能性がチェックされる。
【００８２】
前述のように、ｙモジュール・エージェントが臨界点に近く、したがってモースモードの場合は、最適な最適化の軌跡ｙ＊を検出するときに、前述の標準的なオイラー順方向積分の手続きは役に立たない。一般に、臨界点の近傍では、ｙモジュール・エージェントｉに対応するヘッシアンのブロック対角（block-diagonal）部分行列Ｈｉ，ｉは固有値が負であり、したがって、計算される方向は降下の方向ではない場合もある。しかし、ヘッシアン部分行列Ｈｉ，ｉは、正定値ヘッシアン部分行列
【数７２】
2007509437.tif 000052
を生成し、状態ベクトルのモースモードの前進積分を実行するように変更できる。
【００８３】
状態ベクトルの前進積分へのモースモードのアプローチの派生は、３つの仮定、すなわち（１）目的関数Ｆ（ｙ（τ），ｕ（τ））は、ｕ（τ）に関してほぼ一定である、（２）ヘッシアン行列Ｈ（ｙ（τ），ｕ（τ））は対角中心であり、これはｙモジュール・エージェントが区分化され、ある程度独立していることを示している、（３）状態ベクトルｙｉ（τ）は、ｙモジュール・エージェントｉがモースモードの間はほとんど変化しない、に基づいている。標準的な方法で計算されたρｉ（ｙ（τ），ｕ（τ））の値で決定されるように、
【数７３】
2007509437.tif 000053
で、ｙモジュール・エージェントｉがモースモードすなわち局所的な臨界状態に入ると仮定する。前述の第１と第３の仮定に従って、領域内の点
【数７４】
2007509437.tif 000054
に関して、テーラー近似式（Taylor approximation）により、目的関数Ｆは以下に示すように近似できる。  
【数７５】
2007509437.tif 000055
３次の項を切り捨てると、目的関数は次のように近似できる。  
【数７６】
2007509437.tif 000056
これで、領域内の点
【数７７】
2007509437.tif 000057
を中心とするブロック対角のヘッシアン部分行列Ｈｉ，ｉは、弱いモース変換を使用して次のように正定値行列に変換できる。
【００８４】
まず、図４８に示すように、ヘッシアン行列は次のようにスペクトル分解によって分解できる。  
【数７８】
2007509437.tif 000058
ＮｘＮの行列Ρ  ４８０２は、ヘッシアン行列の固有値を列として備えており、ＮｘＮの対角行列Λ  ４８０４には対応するヘッシアン行列の固有値が含まれる。
【００８５】
ヘッシアンを変換する１つのアプローチは、行列Ρの代わりに新しい行列
【数７９】
2007509437.tif 000059
を作成することである。ただし、この行列の列ベクトルは次のように定義される。  
【数８０】
2007509437.tif 000060
さらに、対角エントリ
【数８１】
2007509437.tif 000061
を用いて、新たな行列
【数８２】
2007509437.tif 000062
が定義される。  
【００８６】
以前の結果を用いると、正の定符号（positive definite）である変更されたヘッシアン行列は、変更したＰ行列および変更したΛ行列から次のようにして得られる。  
【数８３】
2007509437.tif 000063
さらに、変更されたヘッシアン行列の逆行列は、次のようにして計算できる。  
【数８４】
2007509437.tif 000064
逆行列を必要とするのは、変換されたヘッシアンが特異の可能性があるためである。
【００８７】
ヘッシアンを変換する第２のアプローチは、強モース変換（Strong-Morse Transform）として知られており、モースモードの順方向積分ステップで使用される。この方法では、Λを変更した新しい行列が次のようにして生成される。  
【数８５】
2007509437.tif 000065
そして、変更されたヘッシアンは以下のように算出される。  
【数８６】
2007509437.tif 000066
この変換されたヘッシアンは正則（non-singular）であり、したがって、変換されたヘッシアンの逆行列は次のように明確に定義される。  
【数８７】
2007509437.tif 000067
代替の実施形態では、弱モース変換（weak Morse transform）をモースモードでない点、すなわち臨界点から遠い点に使用できる。
【００８８】
変換されたヘッシアン、
【数８８】
2007509437.tif 000068
を使用すると、目的関数の近似は次のようになる。  
【数８９】
2007509437.tif 000069
状態ベクトルｙｉ（τ）の第ｉ区分に関してこの近似目的関数の導関数をとると、次のようにベクトルＧの第ｉ区分が生成される。  
【数９０】
2007509437.tif 000070
さらに、τｋで評価されたベクトルＧの第ｉ区分は次のようになる。  
【数９１】
2007509437.tif 000071
【００８９】
したがって、ヘッシアン部分行列は次のように指定される。  
【数９２】
2007509437.tif 000072
繰り返すが、非対角ヘッシアン部分行列は強モース変換では変更されない。ヘッシアンは対角中心であり、非対角部分行列は目的関数の近似にはほとんど影響しないためである。前述の仮定２を使用すると、τｋにおけるブロック対角ヘッシアン部分行列は、τｉ（上部にハット）におけるブロック対角ヘッシアン部分行列にほぼ等しい。  
【数９３】
2007509437.tif 000073
同様に、第３の仮定は以下を意味している。  
【数９４】
2007509437.tif 000074
したがって、τｋで評価された近似目的関数
【数９５】
2007509437.tif 000075
のベクトルＧは、次のように適切に近似される。  
【数９６】
2007509437.tif 000076
で評価された近似目的関数のヘッシアン部分行列は、次のように近似される。  
【数９７】
2007509437.tif 000077
このような近似によって、エージェントｉに関するｗベクトル区分と局所的なρの近似計算ｗｉとρｉが次のようにして得られる。  
【数９８】
2007509437.tif 000078
【００９０】
図４６のステップ４６３０で、束縛変数の係数を調整してもよい。図１９に示すように、束縛変数の係数は連続的に調整されるが、本発明の少なくとも１つの実施形態において、調整は２進であり、以下の場合に発生する。  
【数９９】
2007509437.tif 000079
ただし、ｙＮはエンベロープ変数ｚである。  
束縛変数の係数が調整される時点では、当初の小さな定数値から特定の大きな定数に変更される。この点で、状態ベクトルも調整する必要があるので、束縛変数の係数を調整しても、状態ベクトルを実現不可能な点に強制的に配置することはない。状態ベクトルは次のように調整される。  
【数１００】
2007509437.tif 000080
ただし、パラメータ「ipm_c_offset2」の値は１に近い。
【００９１】
ｕ＿ａｎｄ＿ｒモジュールは、内部ループの繰り返しごとに制御ベクトルｕの次の値を計算するのに必要である。この計算は、ｙモジュール・エージェントにおけるｙベクトルの計算によく似ている。前述のように、次の制御ベクトルｕを計算するための繰り返しの式は以下のようになる。  
ｕｋ＋ｌ＝ｕｋ−φｋ
【００９２】
静的な問題を定式化するために、ここでは繰り返しごとにハミルトニアンの凸２次近似（convex quadratic approximation）を使用する。ここでは最小値を表す明示的な式を検出できるので、必ずしもパラメータσを使用する必要はない。ｕを計算するには、唯一のステップが必要である。
【００９３】
ここで、ｕｋ＝ｕ|τｋを、関数Ｈ（ｙ，ｕ，ｐ）の点とする。ｙｋとｐｋが固定された場合に、二次テーラー展開を使用したｕｋの周辺の二次近似は、次のようにして計算できる。  
【数１０１】
2007509437.tif 000081
ここで、Ψは∂Ｈ／∂ｕとして定義される。これは２次関数であるが、必ずしも凸ではない。∂Ψ／∂ｕが必ずしも正定値（positive definite）ではないためである。モース変換を使用して、∂Ψ／∂ｕを正定値の凸近似
【数１０２】
2007509437.tif 000082
に変換できる。これは、前述のヘッシアン行列の変換と同様である。
【００９４】
これで、ハミルトニアン行列を次のように近似できる。  
【数１０３】
2007509437.tif 000083
Δｕ＝ｕ−ｕｋと定義すると、次のようになる。  
【数１０４】
2007509437.tif 000084
【００９５】
以下の説明では、ｙｋ、ｕｋ、ｐｋへの依存性については明確に示していない。上の関数の最小値は、目的関数の勾配を０に設定し、Δｕについて解くことによって検出できる。したがって、
【数１０５】
2007509437.tif 000085
となり、次の式が得られる。  
【数１０６】
2007509437.tif 000086
ここで、上付き文字Ｒは前述のように反転を表す。ここで、
【数１０７】
2007509437.tif 000087
と定義すると、次の式が得られる。  
【数１０８】
2007509437.tif 000088
これで、次の制御演算の反復更新が可能になる。  
ｕｋ＋１＝ｕｋ＋Δｕ
しかし、Ψおよび∂Ψ／∂ｕについて解くのは複雑な計算プロセスである。ハミルトニアンは次のように記述される。  
【数１０９】
2007509437.tif 000089
【００９６】
チェインルール（chain rule）と、制御演算ではｐが定数として扱われるという事実とを使用すると、α＝０，．．．，Ｍ−ｌにおいて、
【数１１０】
2007509437.tif 000090
となり、α，β＝０，．．．，Ｍ−１において、
【数１１１】
2007509437.tif 000091
となる。
【００９７】
ここで、α＝０，．．．，Ｎ−１において、次の定義、
【数１１２】
2007509437.tif 000092
により、このｕに関する１次導関数は次のように表される。  
【数１１３】
2007509437.tif 000093
ただし、α＝０，．．．，Ｎ−１、β＝０，．．．，Ｍ−１である。
【００９８】
ヘッシアン行列の逆行列Ｈ−１は次のように定義される。  
Ｈ−１Ｈ＝Ｉ
さらに、行列の要素について次のように記述できる。  
【数１１４】
2007509437.tif 000094
【００９９】
ただし、β＝０，．．．，Ｎ−１であり、さらに
【数１１５】
2007509437.tif 000095
である。両辺で積の法則を使用してｕに関する導関数をとると、次の式が得られる。  
【数１１６】
2007509437.tif 000096
ただし、α，η＝０，．．．，Ｎ−１、γ＝０，．．．，Ｍ−１である。移項すると、次のようになる。  
【数１１７】
2007509437.tif 000097
両辺にＨ−１をかけて次の式を得る。  
【数１１８】
2007509437.tif 000098
式を簡素化すると、次のようになる。  
【数１１９】
2007509437.tif 000099
Ｈ−１の要素を表す明示的な式を次に示す。  
【数１２０】
2007509437.tif 000100
Ｈ−１の式を
【数１２１】
2007509437.tif 000101
に代入すると次の式が得られる。  
【数１２２】
2007509437.tif 000102
ただし、α＝０，．．．，Ｎ−１、β＝０，．．．，Ｍ−１である。
【０１００】
２次導関数は次のようになる。  
【数１２３】
2007509437.tif 000103
【０１０１】
Ｆの勾配とヘッシアン行列はｕに関して線形なので、これらの項のｕに関する２次導関数は０であるという事実によって式が簡素化されることに留意されたい。Ｈ−１の式を上の式に代入すると次の式が得られる。  
【数１２４】
2007509437.tif 000104
【０１０２】
これは次のように簡素化できる。  
【数１２５】
2007509437.tif 000105
【０１０３】
次に、分母（denominator）ρが次のように定義される。  
【数１２６】
2007509437.tif 000106
ｕに関する１次導関数と２次導関数は、積の法則を使用して次のように表される。  
【数１２７】
2007509437.tif 000107
ただし、α＝０，．．．，Ｍ−１、
【数１２８】
2007509437.tif 000108
ただし、α，β＝０，．．．，Ｍ−１。
【０１０４】
次に、形式的な分散アルゴリズムについて説明する。  
【数１２９】
2007509437.tif 000109
の両辺に
【数１３０】
2007509437.tif 000110
をかけて次の式を得る。  
【数１３１】
2007509437.tif 000111
ｒ番目のエージェントの項における関係は次のように表される。  
【数１３２】
2007509437.tif 000112
ただし、ｉ＝０，．．．，Ｍ−１。移項して次の式を得る。  
【数１３３】
2007509437.tif 000113
【０１０５】
以下の表記  
【数１３４】
2007509437.tif 000114
を使用して、両辺に
【数１３５】
2007509437.tif 000115
の逆行列をかけると次の式が得られる。  
【数１３６】
2007509437.tif 000116
したがって、次の式が得られる。  
【数１３７】
2007509437.tif 000117
したがって、制御を次のように更新できる。  
ｕｒｋ＋１＝ｕｒｋ＋Δｕｒ  。
【０１０６】
ここで、エージェントの項に以下に示すハミルトニアンの定義を記述する。  
【数１３８】
2007509437.tif 000118
これで次の式が得られる。  
【数１３９】
2007509437.tif 000119
【数１４０】
2007509437.tif 000120
ここで、
【数１４１】
2007509437.tif 000121
（ただし、α＝０，．．．，ｎｉ）、および
【数１４２】
2007509437.tif 000122
を考慮すると、エージェントＩのｗベクトルの各要素は次のように表すことができる。  
【数１４３】
2007509437.tif 000123
【０１０７】
前の式にＨ−１の定義を代入すると、次のようになる。  
【数１４４】
2007509437.tif 000124
移項すると次の式が得られる。  
【数１４５】
2007509437.tif 000125
ＧとＨのｕに関する２次導関数は０なので、次のように表される。  
【数１４６】
2007509437.tif 000126
【０１０８】
ただし、
【数１４７】
2007509437.tif 000127
ここで、項を分割し、終わりの２項の添字を変更すると、次のようになる。  
【数１４８】
2007509437.tif 000128
したがって、次の式が得られる。  
【数１４９】
2007509437.tif 000129
【０１０９】
さらに、ＧとＨのｕに関する２次導関数は０であるという事実を使用して、次の式を得る。
【数１５０】
2007509437.tif 000130
最終的に次の式が得られる。  
【数１５１】
2007509437.tif 000131
さらに、この導関数は次のようになる。  
【数１５２】
2007509437.tif 000132
【０１１０】
共役変数の伝搬は、状態変数と制御変数の伝搬とは異なっている。状態変数と制御変数は、密結合し、共に前方に伝搬する。対照的に、共役変数は状態変数と制御変数の結合とは関係なく、適切な時期に後方に伝搬する。このためには、ｙモジュールとｕ＿ａｎｄ＿ｒモジュールで計算した値のデータベースが必要である。ｐの微分方程式はｐに関して線形なので、解は積分でなく補間法によって近似できる。多くのアプリケーションで、ベクトルｐの値は収束限界周辺の過渡現象（transient around the convergence horizon）を除いてほとんど一定である。ｐの式における行列Ａもｐに依存しない。
【０１１１】
ｐを計算するために、ｙモジュールとｕ＿ａｎｄ＿ｒモジュールで計算されたｙとｕの軌跡を表す一連のノードの値が作成され、これをｐモジュールで使用する。まず終点の条件ｐ（Ｔ）が計算され、次にｙの初期の計算まで後方に補間するループの中で、ｐの初期の値が計算される。  
【数１５３】
2007509437.tif 000133
【０１１２】
実用的にするために、ここでは行列Ａがほとんど一定であると仮定される場合は、τｋ−１−τｋ（ただし、τｋ−１＜τｋ）のサイズの増分に問題を分割する。次に、後退積分を実行すると、τｋが指定された場合にτｋ−１における解は次のように表される。  
ｐ(τｋ−１)＝exp(Ａ(τｋ−１−τｋ))ｐ(τｋ)
【０１１３】
この式の指数は、以下の近似によって解決される。  
【数１５４】
2007509437.tif 000134
ただし、Ｎは近似の次数である。
【０１１４】
修正（Repair）
ここまで、一般的な最適化の特定の計算方法について説明してきた。しかし、現実世界の最適化における付加的な側面は、実際の時間区間にわたって繰り返される最適化に関連する。たとえば、メーカーなどの企業では、最適化の技術を使用して、従業員のスケジューリング、パーツや原料の注文、パ―ツ、原料、中間製品の製造工場内での配布、そして特に顧客への製品の配送を含む製造工場の運転に関するさまざまな側面をすべて最適化しようとする。一般に、こうした企業の最適化を行うための非常に高次元の問題領域では、製造工場運転中に発生する管理されない多くのイベントや変化によって、時間が経つと共に変化する条件に合わせて製造工場の制御を再調整するために、最適化を何度も実行しなければならない場合がある。最適化以外に、同様の問題は他の多くの計算タスクにも発生する。たとえば、複雑な系のモデル化では、系に関してランダムに発生する様々なイベントや混乱に対応するために、系の状態に関する中間的な計算を導入し、変化する環境への系の実際の時間応答をそのモデルが反映するようにしなければならない場合がある。
【０１１５】
図４９は、こうした反復最適化の戦略を示している。図４９と後続の図５０および５２〜５６において、横軸４９０２は時間に対応する。図４９では、状態変数と制御変数を含む系の準最適な状態関数が、図４９に縦の破線で示す各時点（たとえば縦の破線４９０４）で、最適化および／またはモデル化の技術を使用して計算される。計算された状態関数は、図４９には曲線の一部分として表されている。たとえば、曲線の一部分４９０６は、時刻４９０８の点で開始されており、この時点で状態関数の計算が開始され、別の時刻４９１０の点まで展開されている。もちろん、計算された状態関数は、現実の世界に適用する上で、超次元空間における複雑な制御の軌跡であるが、説明を明確にするために、曲線の一部分４９０６のような２次元の制御の軌跡を使用している。斜めの破線（たとえば斜めの破線４９１２）は、前述の連続型反復変数τに関する積分を使用した状態関数の計算を表しており、状態関数を計算する最適化技術の計算限界Ｔに関連する特定の長さの時間４９１４を必要とする。新しい状態関数、たとえば、時刻４９１８で開始され、時刻４９２０で終了する時間の範囲に適用できる新しい状態関数４９１６が計算されると、この新しい状態関数は直前の状態関数（この例では状態関数４９０６）と一部の重複する部分で比較的接近してはいるが、一致してはいない。このように、連続して計算された２つの状態関数の間には、ずれ４９２２または不連続が存在する。この不連続は、前の状態変数が計算された時刻と現在の状態変数が計算された時刻との間に発生した系の条件や環境の変化によって発生する。この場合は、時刻４９０８で開始され、時刻４９１８で終了する時間区間にイベントや変化が発生している。
【０１１６】
残念ながらこうした不連続は系の混乱を表している。たとえば、製造工場において、不連続は従業員のスケジュール変更、配送用トラックの出発時刻と行き先の変更、およびその他の同様の変更を表していてもよい。しかし、複雑な系は慣性（inertia）をもつ傾向があるので、大きな混乱を克服するのは困難であり、前回計算した最適な状態関数が最適でなくなる段階まで、自らが系の環境をある程度変更し、複雑なシステムをさらに混乱させるリプルやカスケードの効果を発生するという問題を克服するには、長い時間と多くの資源が必要になる場合がある。
【０１１７】
最適化の繰り返しによって発生する不連続を改善する１つの可能なアプローチは、最適化を実行する間隔を短くすることである。図５０は、図４９に示す状態関数計算による最適化の間隔を短縮したものを示す図である。状態変数をより効率的に計算することによって、図５０では状態関数の計算を表す破線の傾きが大きくなっており（たとえば破線５００２）、最適化を短い間隔ΔＴ  ５００４で実行してもよい。しかし、その結果として混乱の頻度が増し、混乱の頻度が増したこと自体によってリプルやカスケードの効果が発生し、場合によっては不安定性がますます増大することにもなり得る。一部の事例では、非常に短い時間区間に特定の種類のイベントが同時に発生したことによる不足減衰（under damping）の結果として、より大きなずれが発生することさえある（たとえば、ずれ５００６）。
【０１１８】
図５１は、最適化の繰り返しによって発生するもう１つの結果の類推を示している。図４９〜５０では、新しい状態関数は、前に計算された状態関数による制御を終了し、新しく計算する状態関数による系の制御を開始する点を除いて、前に計算された状態関数に関係なく新たに計算されることを想定していることに留意されたい。しかし、前の状態関数計算の履歴を考慮しないと、長い時間が経過した場合に、極めて非最適な状態関数の解が得られる可能性がある。同様に、図５１に示す単純な地形図で表される地域を点Ａ  ５１０２から点Ｂ  ５１０４まで移動する問題について考察する。点Ａから点Ｂに至る第１のパス５１０６は、ある人がコンパスだけを頼りにたどるであろうパスを表しており、コンパスの針の動きが止まるのを待たずに、瞬間的なコンパスの読みだけに基づいて頻繁に軌道修正を繰り返している。こうしたパスは、点ＡとＢの間の移動に関して、時間的にも距離的にも全くの非最適解を表している。
【０１１９】
第２のパス５１０８は、コンパスの動きが落ち着いてから読む必要があることを認識している人、がけに磁鉄鉱が露出している場合などの局所的な条件によって混乱が発生する可能性があることを認識している人、および障害物が現れた場合は、コンパスの示す方向でなくその障害物に慎重に対処することを認識している人がたどるであろうパスを表している。第２のパス５１０８は、点ＡとＢの間の移動に関して、はるかに効率的かつ効果的なパスを表しているのは明らかである。
【０１２０】
第３のパス５１１０は、地形数学（topographic math）を理解でき、コンパスを正しく使用でき、川にかかった橋は目的地への重要な中間点であることを認識している人、森の道に沿って歩く方が道を開拓するよりはるかに容易なことを認識している人がたどるパスである。図５１から類推すると、過去の情報に基づくことができる人、単純で限られた現在の状態以外の情報に基づくことができる人、追加の情報に基づいて将来的な解を予測できる人は、点ＡとＢの間の最適なパスを検出する可能性の高い人である。図４９〜５０に示す反復最適化問題に戻り、局所的、瞬間的な状態関数を考慮するのみでなく、前に計算された状態関数を使用してより長い期間にわたってより多くの情報を取得する方法を検出することが望ましい。
【０１２１】
図４９〜５０に関連する前述の考慮事項の結果として、長期にわたる履歴情報と最新の予測を使用してある程度連続型最適化のアプローチは、最適化を目的とする本発明の実施形態に追加される機能、「修正に基づく連続型最適化のアプローチ」（ＲＣＯ：Repair-based Continuous Optimization）を表している。図５２は、変動時間枠（sliding window）の概念を示している。図５２に示すように、時刻ｔ１−ΔＴ  ５２０４で始まり、時刻ｔ１＋Ｔ−ΔＴ  ５２０６で終了する時間長Ｔの前の時間枠５２０２は、状態関数を計算するための参照フレームとして使用される。この時点から、時間長ΔＴ  ５２０８が経過している。したがって、時間長Ｔの新しい現在の時間枠５２１０は、現在の時刻ｔ１  ５２１２における参照フレームを提供する。
【０１２２】
図５３に示すように、前の時間枠５３０２で計算された状態関数は時刻ｔ１−ΔＴ  ５２０４からＴ１−ΔＴ＋Ｔ  ５２０６に展開される。図５４は、前に計算された状態関数の現在の時間枠への展開を示している。図５４に示すように、前に計算された状態関数５３０３は、５４０２から現在の時間枠の終わり（時刻ｔ１＋ΔＴ  ５４０４）まで展開される。図５５は、現在時点における考慮事項、過去のイベントに関する考慮事項、および現時点で前に計算された状態関数が陳腐化し、最適でなくなる可能性のある環境の変化、および現在のウィンドウで表される現在の参照フレームで前の状態関数の情報を新しい最適な状態関数の検出に組み込む方法を示している。
【０１２３】
図５５において、縦の線５５０２〜５５０５は、状態関数の計算時間中、すなわちｔ１−ΔＴおよび現在時刻ｔ１に発生したイベントを表している。こうしたイベントによって、系の状態や環境が変化し、前に計算された状態関数が最適でなくなる可能性がある。たとえば、製造工場で、ストライキやインフルエンザ・ウィルスによって多くの従業員が製造工場内の作業を実行するようにスケジューリングできなくなる可能性がある。したがって、前に計算された最適なスケジューリングは、従業員不足に直面した場合は完全な次善最適解であろう。図５５に示すように、状態関数の計算の基準となる時間依存パラメータに関するさまざまな予測は、時刻ｔ１では変化している可能性がある。こうした予測の変化は、破線５５０６で表されている。たとえば、製造工場のコンテクストにおいて、主要な顧客がこの顧客が通常購入する数の２倍の製品を緊急に要求してきたり、競合他社で製造された新しい製品の出現によってこの製造工場で製造する製品の需要が突然減少したりする可能性がある。このように、前に予測した製品需要が現在および将来の需要を反映しない場合がある。
【０１２４】
さらに、状態関数の基準となるさまざまな考慮事項やパラメータに対して当初指定された優先順位が変化する可能性がある。図５５の点線５５０８は、現在の参照フレームで発生した優先順位の変化を前の参照フレームに関連して示している。たとえば、木曜日に状態関数が計算された時点で、時刻ｔ１−ΔＴからｔ１までの間に、差し迫った道路の閉鎖に関する情報が発表されていなかったとする。木曜日の道路の閉鎖に関する情報によって、南側に住む従業員が通勤時に橋の閉鎖の影響をほとんど受けない時間帯で作業できるスケジューリングの優先順位が大幅に向上する。
【０１２５】
さらに、系の動作における重大な混乱を回避するために、系は新しく計算された状態関数が特定の区間にわたって同じ状態を示しているかどうかにかかわらず、直前に計算された系の状態に従う必要がある。たとえば図５５で、それぞれＴ２  ５５１４とＴ４  ５５１６で始まり、それぞれＴ３  ５５１８とＴ５  ５５２０で終了する２つの時間枠５５１０と５５１２は、新しく計算された状態関数内に状態関数の古い値５５２２と５５２４を保持する必要がある。このように、前の状態関数の計算で得られる履歴の情報は、現在の参照フレームに保持される。このようなフェンス区間は、前に計算された状態関数に基づいて確立するシステム制御における大きな障害や混乱を避けるために必要であろう。
【０１２６】
イベントに関する情報、環境や優先順位の変化に関する情報、および前に計算された特定の状態を保持する必要があるかどうかを使用して、新しい状態関数を計算せずに、前に計算され、展開された状態関数を変更し、現在の時間枠に関する準最適解を生成することができる。図５６は、現在の時間枠に関する準最適な状態関数を得るための、前に計算された状態関数の修正を示す図である。図５６に示すように、前に計算された状態関数の修正には、現在の時間枠内の時刻における、前に計算された状態関数と現在の時間枠に関する新しい準最適な修正された状態関数との差分の計算が含まれる（たとえば差分５６０２）。差分５６０２のような差分の計算では、前の時間区間ｔ１−ΔＴ  ５２０８に発生したイベント、環境や優先順位の変化（図５５の５５０８と５５０６）、特定の時間区間５５１０と５５１２に関して前に計算された状態５５２２と５５２４を保持する必要があるかどうかが考慮される。したがって、修正に基づく連続型最適化技術によって、系の制御における混乱や障害を最小化し、不連続な制御を回避するように最適化全体を繰り返し実行することなく、ある程度連続型状態関数の再計算が可能になる。
【０１２７】
以下は、本発明の１つの実施形態を表す修正方法の、より形式的な記述である。制御最適化問題は、その解が前述の修正方法で修正される制御状態関数を提供し、以下のように指定される。  
【数１５５】
2007509437.tif 000135
【０１２８】
ただし、制約
【数１５６】
2007509437.tif 000136
を受け、ここで、ＶはＲｒのコンパクト・サブセット（compact subset）であり、ｄ（ｔ）は駆動（driving）関数（たとえば、決定論（deterministic）または予測需要（forecasted demand））である。  
【数１５７】
2007509437.tif 000137
および｛ｇｊ，ｊ＝１，．．．，ｍ｝は、十分になだらかである。以下のようなスカラー変数ｘｎ（ｔ）≧０が導入される。  
【数１５８】
2007509437.tif 000138
これで、最適制御問題は次のようになる。  
【数１５９】
2007509437.tif 000139
ただし、ｔ∈［ｔ１，ｔ１＋Ｔ］の場合
【数１６０】
2007509437.tif 000140
を受け、ここで、
【数１６１】
2007509437.tif 000141
である。こうした問題の解ｖ＊（ｔ）は、「最適制御（optimal control）」と呼ばれており、対応する軌跡ｘ＊（ｔ）は最適な軌跡（optimal trajectory）と呼ばれている。前述の制御状態関数（たとえば、図４９の４９０６や４９１６）は、ｘ（ｔ）とｖ（ｔ）（ｔ∈［ｔ１，Ｔ］）との組合せである。
【０１２９】
あらかじめ計算された制御状態関数は、図５２〜５６で使用して表すと、次のｖ（ｔ１−ΔＴ）＊（ｔ）およびｘ（ｔ１−ΔＴ）＊（ｔ）の組合せとなる。計算された差分（たとえば、図５６の５６０２）は、「修正（the repair）」と呼ばれ、δｘｔ１ｌ（ｔ）およびδｖｔ１ｋ（ｔ）の組合せとして定義される。フェンス・インディケータ（Fencing indicator）関数は、制御状態変数のそれぞれについて、ｔが修正された制御状態関数に偏差が許容されない区間内にあるかどうかを指定する。  
【数１６２】
2007509437.tif 000142
【０１３０】
時間区間［ｔ１，ｔ１＋Ｔ］内の有限個のフェンス区間（fenced intervals）を想定すると、（ｉ＝１，．．．，ｎ）（ただし、ｎはモデルの次元）はモデルの次元であり、増大する状態の軌跡δｘｔ１ｉ（ｔ）は次のように指定される。  
【数１６３】
2007509437.tif 000143
【０１３１】
ここで、（ｔ１−Δｔ）＊の表記は、偏導関数が参照軌跡に沿って評価されることを示している。
【０１３２】
集合｛[ｔｓｉ，ａ、ｔｓｉ，ｂ]、ｓ＝１,...,Ｍｉ｝は、状態ｉ＝１,...,ｎのフェンスのある時間の部分区間である。時間区間［ｔ１−ΔＴ，ｔ１］内に発生するすべてのイベントＡｅｉδ(ｔ−ｔｅ)は、次のようにｔ１に発生する唯一のイベントＡｔ１ｉδ(ｔ−ｔ１)に近似される。  
【数１６４】
2007509437.tif 000144
δｘｔ１ｉ（ｔ１）＝ｘｔ１ｉ（ｔ１）−ｘ（ｔ１−ΔＴ）＊(ｔ)−Ａｔ１ｉは、δｘｔ１ｉ（ｔ１）の初期条件を指定する。ただし、ｖ（ｔ１−ΔＴ）＊（ｔ）およびｘ（ｔ１−ΔＴ）＊（ｔ）は、［ｔ１−ΔＴ，ｔ１＋Ｔ−ΔＴ］上で定義される。前述のように、制御状態関数は区間［ｔ１，ｔ１＋Ｔ］に展開する必要がある。展開は一定の制御によると仮定し、また展開内のフェンスは許可されないと仮定すると、次式のようになる。  
【数１６５】
2007509437.tif 000145
ここで、ｔ∈［ｔ１＋Ｔ−ΔＴ，ｔ１＋Ｔ］（ｉ＝１,…,ｎ）であり、
【数１６６】
2007509437.tif 000146
である。
【０１３３】
区間［ｔ１，ｔ１＋Ｔ］に関する動力学（dynamics）は、次のように与えられる。  
【数１６７】
2007509437.tif 000147
Φ（ｘ（ｔ１＋Ｔ））で
【数１６８】
2007509437.tif 000148
を表すと、２次までの展開によって次の式が得られる。  
【数１６９】
2007509437.tif 000149
【０１３４】
修正の目的は、上の基準を最小化しながら、同時にノミナル（nominal、名目上の、）軌跡に関する変化を最小化することである。この基準の組合せは次のとおりである。  
【数１７０】
2007509437.tif 000150
ここで、ＱとＲは企業の基準（enterprise criterion）を満たすことと変化を最小化することとの兼ね合いを定義するためにユーザーが指定する定数の正定値行列である。
【０１３５】
修正問題の制約は次のとおりである。  
【数１７１】
2007509437.tif 000151
この式を１次まで展開すると次のようになる。  
【数１７２】
2007509437.tif 000152
【０１３６】
制約の不等式を満たすために、次の式が与えられる。  
【数１７３】
2007509437.tif 000153
【０１３７】
展開された制約の第２項と第３項は次を満たす必要がある。  
【数１７４】
2007509437.tif 000154
【０１３８】
メタ制御と階層的プロセスの構築
これまでの記述は、最適化と後続の点での最適化法によって生成された制御状態関数の調整に焦点を合わせている。しかし、本発明には、複雑で階層的な計算システムとモデルを構築するため、さらに簡単か複雑かにかかわらず一般的な計算を制御するための、より一般的で強力な方法が含まれる。前述の計算制御によるこの一般的な方法の例を示し、当初の最適化問題に追加の制約を加えてメタレベルの最短時間制御問題の目的関数Ｊを生成し、さらにメタレベルの問題に関するハミルトニアンを最小化することによって準最適制御ポリシーを生成する。
【０１３９】
図５７は、計算プロセスのメタ制御の１つのレベル示している。図５７では、最初の計算問題または計算モデルである５７０２は、メタレベル軸５７０４の最も低いレベルに示されている。計算問題、またはモデルは離散的であり、反復
【数１７５】
2007509437.tif 000155
によって形式的に説明できる。前述のように、離散的な計算問題またはモデルに連続型反復変数のパラメータ化を導入することによって、当初の計算問題を連続的とみなすことができ、問題を解決するために強力で連続型計算方法を適用できる。したがって、図５７に示すように、この問題は面５７０６では離散的に思われるが、別の面５７０８では連続型（continuous）ｄｘ（τ）／ｄ（τ）＝ｆ（ｘ（τ）、ｕ（τ））とみなされる。離散から連続への変換は、メタ制御の方法を有効にする重要な機能である。問題の連続型ビューを使用することにより、追加の制約を導入でき、前述の最適化問題で構築された最短時間制御問題などの新しい問題を作成できる。この制約を追加する機能の重要な側面は、新しく追加された制約が限定的ではないこと、または当初の問題への関連付けが必要ないことである。また、前述の実施形態において、初期の問題には最適化問題に対する最適制御解の検出が含まれるが、追加の制約は当初の最適化問題において時間的に効率のよい計算に向けて計算を操作することであり、準最適制御ポリシーの検出に関連する。その他の種類の制約を追加してもよい。たとえば、最適化を実行する曜日を操作するメタレベルの制約を追加してもよい、時間とメモリの両方を効率的に使用して計算するための制約を追加してもよい。このように、追加の制約は当初の問題に重ね合わせるメタレベルのビューを表している。新しい制約を追加できることによって、たとえば、周知のトップダウン型アプローチで問題を解決するといった複雑なソフトウェアの設計で使用される概念的な設計のレベルと同様に、階層的な計算による問題解決方法を容易に適用できる。
【０１４０】
図５７では、大きな矢印５７１０で表される変換の最初の部分として、制約の追加が示されている。当初の問題全体が、新しい問題における制約の内部と新しい目的関数の項内の両方で、新しい問題に伝えられていることに留意されたい。したがって、新しい問題は当初の問題を含んでおり、当初の問題の解は新しい問題の解の部分集合になっている。したがって、メタレベルではあるが、新しい制約と当初の問題とは全く異なるビューに焦点を合わせることができ、当初の問題は新しい問題として符号化される。変換の第２の要素では、ハミルトニアンが順に使用され、メタ制御ポリシーを計算するための離散型（discrete）の計算方法を決定するための式を提供する。離散型のメタ制御ポリシーの式５７１２は、新しい連続する繰り返し変数σを導入することで、それ自体を連続型の問題５７１４とみなすことができる。
【０１４１】
図５８は、複数レベルのメタ制御および／または複雑な計算プロセスの階層構造を示している。図５８に示すように、当初の離散問題５８０２からのメタレベルの変換を繰り返すことで、メタレベルの軸５８１０に沿って任意の数のメタレベル問題が得られる。前述のように、それぞれの新しいメタレベルにより、当初の問題および低レベルのメタレベルを高レベルのメタレベルへと伝搬させながら、新しい制約を追加し、新しい問題を生成できる。例として、現在科学界で大きな関心がよせられている、生きた細胞をモデル化するという計算によるモデル化の問題を考える。細胞は数多くの概念的なレベルとして階層的に見ることができる。すなわち、量子力学および統計力学によって支配される分子レベルから始めて、超分子複合体（super  molecular complexes）、巨大分子構造（macromolecular structures）、およびコンパートメント（compartments）、細胞小器官（organelles）、遺伝子制御（genetic control）、環境の影響、および他の無数の概念レベルへと向上し、それぞれは異なる観点でとらえることができ、さまざまな力と原理で操作される。しかし、各レベルは細胞全体に寄与し、分子レベルのイベントの作用や状態は上記の概念レベルに深く影響している。たとえば、特定の染色体内における核酸塩基から別の核酸塩基への自然変換は、遺伝子の変化を引き起こし、次に重大な組織構造または一時的な制御構造の構成要素である遺伝子産物（gene product）の変化を引き起こす。これは、さらに細胞全体の存続に重大な影響を与える可能性がある。さらに、各レベルの問題のさまざまな部分は、本質的に整数、ブール値、または他の型の変数で表現できるが、問題の他の部分は連続変数の項で定義できる。
【０１４２】
本発明の方法は、生きた細胞やリアルタイムの分散オペレーティング・システムなどのシステムの複雑な計算モデルを開発するのに適用できる。前述の非実数変数を実数変数に変換し、さらに浮動小数点変数に変換する方法は、さまざまな下位の問題の本質的な表現を、本発明の方法を適用できるすべて浮動小数点からなる問題空間に組み込む方法を提供する。さまざまなメタレベルのビューを定義できることにより、低レベルのビューを伝搬させながら、高レベルのメタレベルで低レベルのビューを明示的に指定する必要なしに、細胞のモデル化の問題から低レベルのビューの複雑性をすべて保持する階層的なモデルの作成が可能になる。
【０１４３】
より形式的には、次のように反復アルゴリズムが提供される。
【数１７６】
2007509437.tif 000156
ここで、ｕｋは制御変数、ｘｋは決定変数である。このアルゴリズムは次の式に収束する。  
【数１７７】
2007509437.tif 000157
【０１４４】
このアルゴリズムの詳細は、関数ｆ(チルダ)とパラメータｕｋでコード化される。
【０１４５】
連続化（continualization）の手続きでは、反復はｘ（τ）とｕ（τ）に関する
【数１７８】
2007509437.tif 000158
となるような制御された微分方程式に変換される。変数τは、離散型の反復変数ｋに対応する連続型反復変数である。この微分方程式は、すべてのｕ（□）∈Ｕ（□）について次の形をとる。  
ｄｘ（τ）／ｄτ＝ｆ（ｘ（τ），ｕ（τ））
ここで、Ｕ（□）は実現可能なパラメータ値の集合である。
【０１４６】
最適制御問題は次のように表現される。  
【数１７９】
2007509437.tif 000159
【０１４７】
関数φとΨは、アルゴリズム設計者がメタ目的（meta objectives）を達成するために選択する。
【０１４８】
ポントリャーギン（Pontryagin）の最小原理と、
【数１８０】
2007509437.tif 000160
で与えられるハミルトニアン（ここで、ｐ（τ）は共役（costate）である）を使用すると、ｘ＊（τ），ｕ＊（τ），∈［０，Ｔ）が最適であるために必要な条件は次のとおりである。
【数１８１】
2007509437.tif 000161
ここで、ハミルトニアンは次のようにｕ＊（τ）におけるｕの関数としての絶対的な最小値をとる。  
Ｈ(ｘ＊(τ)，ｕ＊(τ)，ｐ(τ))≦Ｈ(ｘ＊(τ)，ｕ(τ)，ｐ(τ))，
∀ｕ＊(τ)∈Ｕ(τ)，τ∈［０，Ｔ)。
【０１４９】
収束数列
【数１８２】
2007509437.tif 000162
が、次式により生成される。  
【数１８３】
2007509437.tif 000163
ここで、ｘ（τ）＝ｘ＊（τ）、τ∈［０，Ｔ)のとき、
【数１８４】
2007509437.tif 000164
となる。
【０１５０】
前述の収束数列の生成式は、次に示す当初の反復アルゴリズムに類似している。  
【数１８５】
2007509437.tif 000165
ここで、
【数１８６】
2007509437.tif 000166
は、
【数１８７】
2007509437.tif 000167
に類似している。  
【０１５１】
ｓの連続型パラメータ化をσとして前述の連続化手続きを使用すると、次の式が得られる。    
【数１８８】
2007509437.tif 000168
ここで、ｖ（τ，σ）は
【数１８９】
2007509437.tif 000169
を連続化した形であり、Ｗ（ｘ（τ），ｖ（τ），ｐ（τ））は、
【数１９０】
2007509437.tif 000170
を連続化した形である。  
要約すると次のようになる。
【数１９１】
2007509437.tif 000171
【０１５２】
本発明について、特定の実施形態に関連して説明してきたが、本発明をこの実施形態に限定する意図はない。本発明の精神を逸脱しない変更は、当業者には言うまでもない。たとえば、ほとんど制限なく多種多様である。
【０１５３】
以下の説明において、本発明の理解を深めるために、説明を目的として特定の専門用語を使用した。ただし、特定の詳細が本発明を実施するために必要でないことは、当業者には言うまでもない。他の例では、基盤となる発明からの不要な逸脱を回避するために、周知の回路とデバイスがブロック図の形で示されている。したがって、本発明の特定の実施形態に関する以上の説明は、例示および説明を目的としており、すべてを記述するものではなく、本発明を開示した形式のみに限定するものでもない。以上の教示に関連して、さまざまな変更や変形が可能なことは明らかである。本実施態様は、本発明の原理の最適な説明を提供するために選択され、説明されており、他の当業者はこれを実際に適用することによって、本発明とさまざまな変形を伴うさまざまな実施形態を、予定される特定の用途に適した形で最適に利用できる。本発明の範囲は、前述の請求項とその均等物によって定義されるものである。
【図面の簡単な説明】
【０１５４】
【図ｌＡ】非常に単純な２次元の最適化問題を示す図である。
【図ｌＢ】非常に単純な２次元の最適化問題を示す図である。
【図２】本発明の１つの実施形態を表す最適化の方法およびシステムの最高レベルのビューを示す制御流れ図である。
【図３】モデル作成と前処理を説明するための仮説的な問題を示す図である。
【図４】仮説的な問題のさらなる側面と特性を、図３の描画法を使用して示す図である。
【図５】工場から卸売店、および卸売店から顧客への商品の出荷に関連するコストを示す図である。
【図６】仮説的な問題の付加的な特性を示す図である。
【図７】２か所の製造工場、３か所の卸売店の候補、５か所の顧客による、コスト、容量、需要の数値に関連する単純な仮説的な問題を示す図である。
【図８】図２に示す高レベルの制御流れ図の前処理のステップを示すより詳細な流れ図である。
【図９】図７に示す特定の最適化問題を最初に数理的に表現した記号形式のモデルを示す図である。
【図１０】仮説的な数理モデルから浮動小数点変数のみを使用したモデルへの変換を示す図である。
【図１１】標準の形式に変換した仮説的な問題の数理モデルを示す図である。
【図１２Ａ】変換された等式制約に対する束縛変数の効果を示す図である。
【図１２Ｂ】変換された等式制約に対する束縛変数の効果を示す図である。
【図１２Ｃ】変換された等式制約に対する束縛変数の効果を示す図である。
【図１２Ｄ】変換された等式制約に対する束縛変数の効果を示す図である。
【図１３】等式制約を不等式制約に変換した後の仮説的な問題を記号形式で数理的に表現したモデルを示す図である。
【図１４】束縛変数の効果を図式的に示す図である。
【図１５】エンベロープ変数ｚを導入した後の仮説的な問題の数理モデルを示す図である。
【図１６Ａ】バリア関数の概念を示す図である。
【図１６Ｂ】バリア関数の概念を示す図である。
【図１６Ｃ】バリア関数の概念を示す図である。
【図１６Ｄ】バリア関数の概念を示す図である。
【図１７】バリアの制約を追加した後の仮説的な問題の記号形式による数理モデルを示す図である。
【図１８】初期点（初期点）から局所的最小値に至る最適化の軌跡の方向を示す問題領域内の勾配フィールドの利用を示す図である。
【図１９】最適化の軌跡を示す図である。
【図２０】準最適状態ベクトルの軌跡ｙ＊と準最適制御ポリシーｕ＊およびｒ＊に基づく反復最適化法を図式的に示す図である。
【図２１】反復最適化法全体を図式的に示す付加的な図である。
【図２２】反復最適化のステップ（図２の２０３）を説明する制御流れ図である。
【図２３】椀形の問題領域と多くの中央の突起の面境界を示す図である。
【図２４】図２３と同じ問題領域を示し、臨界点を明示的に示す図である。
【図２５】臨界点の拡大表示を示す図である。
【図２６】第１の検索技術を示す図である。
【図２７】問題領域の第２の検索戦略を示す図である。
【図２８】最適化スレッドのモジュールとモジュール編成を示す図である。
【図２９】スーパーモジュールで保守するデータベースの内容を部分的に示す図である。
【図３０】図２８の描画法を使用して、最適化法の１つの外部ループの繰り返し中に階層構造をなす最適化軌跡のモジュール間で情報が流れる様子を示す図である。
【図３１】図２８の描画法を使用して、最適化法の１つの外部ループの繰り返し中に階層構造をなす最適化軌跡のモジュール間で情報が流れる様子を示す図である。
【図３２】内部ループの次のステップ（図２２の２２０５）におけるデータの流れを示す図である。
【図３３】図３１と３２に示すデータ交換のステップが内部ループの繰り返しによって正常に繰り返される様子を示す図である。
【図３４】図２２に示すステップ２２０９のｐモジュール２８１０で実行する共役ベクトルの計算に含まれるデータ交換を示す図である。
【図３５】ｗベクトルの計算を図式的に示す図である。
【図３６】ヘッシアン行列とｗベクトルの乗算によって得られるＧベクトルを分割された形で示す図である。
【図３７】ヘッシアン行列とｗベクトルの乗算を分解して示す図である。
【図３８】ｙモジュール・エージェントによるヘッシアン行列の逆行列の分解を容易にするＧベクトル区分の計算の数理的な操作を示す図である。
【図３９】第２のｙモジュール・エージェントによるｗベクトル区分ｗ１の計算を示す図である。
【図４０】最高の概念的レベルのスーパーモジュールを示す制御流れ図である。
【図４１】図４０のステップ４００４で呼び出すルーチン「初期点」を示す制御流れ図である。
【図４２】ルーチン「次のＴを選択」を示す制御流れ図である。
【図４３】図４０のステップ４０１２で呼び出すルーチン「評価」を示す制御流れ図である。
【図４４】ルーチン「報告」を示す制御流れ図である。
【図４５】図４０のステップ４０１０で呼び出すルーチン「ｕおよびｙをＴまで積分」を示す制御流れ図である。
【図４６】ルーチン「ｙを積分」を示す制御流れ図である。
【図４７】図４６のステップ４６２５で呼び出すルーチン「ΔＴを計算」を示す制御流れ図である。
【図４８】ヘッシアン行列がスペクトル分解によって分解できることを示す図である。
【図４９】制御問題の最適化法による連続型状態関数の計算を示す図である。
【図５０】図４９に示す状態関数の短縮された区間による計算を示す図である。
【図５１】繰り返しによる現在の状態ベースの最適化から得られる結果の類推を示す図である。
【図５２】変動時間枠の概念を示す図である。
【図５３】前の時間枠で計算されたＴ１−ΔＴからＴＴ−ΔＴまでの状態関数を示す図である。
【図５４】現在の時間枠の時間区間において前に計算された状態関数の展開を示す図である。
【図５５】現在時点における考慮事項と、過去のイベントに関する考慮事項と、現時点で前に計算された状態関数が陳腐化し、最適でなくなる可能性のある環境の変化、および現在の時間枠で表される現在の参照フレームで前の状態関数の情報を新しい最適な状態関数の検出に組み込む方法を示す図である。
【図５６】現在の時間枠の準最適状態関数を得るための前に計算された状態関数の修正を示す図である。
【図５７】唯一のレベルによる計算プロセスのメタ制御を示す図である。
【図５８】複雑な計算プロセスの複数のレベルによるメタ制御および／または階層構造を示す図である。
【図１Ａ】
2007509437.tif 000193
【図１Ｂ】
2007509437.tif 000194
【図２】
2007509437.tif 000195
【図３】
2007509437.tif 000196
【図４】
2007509437.tif 000197
【図５】
2007509437.tif 000198
【図６】
2007509437.tif 000199
【図７】
2007509437.tif 000200
【図８】
2007509437.tif 000201
【図９】
2007509437.tif 000202
【図１０】
2007509437.tif 000203
【図１１】
2007509437.tif 000204
【図１２Ａ】
2007509437.tif 000205
【図１２Ｂ】
2007509437.tif 000206
【図１２Ｃ】
2007509437.tif 000207
【図１２Ｄ】
2007509437.tif 000208
【図１３】
2007509437.tif 000209
【図１４】
2007509437.tif 000210
【図１５】
2007509437.tif 000211
【図１６Ａ】
2007509437.tif 000212
【図１６Ｂ】
2007509437.tif 000213
【図１６Ｃ】
2007509437.tif 000214
【図１６Ｄ】
2007509437.tif 000215
【図１７】
2007509437.tif 000216
【図１８】
2007509437.tif 000217
【図１９】
2007509437.tif 000218
【図２０】
2007509437.tif 000219
【図２１】
2007509437.tif 000220
【図２２】
2007509437.tif 000221
【図２３】
2007509437.tif 000222
【図２４】
2007509437.tif 000223
【図２５】
2007509437.tif 000224
【図２６】
2007509437.tif 000225
【図２７】
2007509437.tif 000226
【図２８】
2007509437.tif 000227
【図２９】
2007509437.tif 000228
【図３０】
2007509437.tif 000229
【図３１】
2007509437.tif 000230
【図３２】
2007509437.tif 000231
【図３３】
2007509437.tif 000232
【図３４】
2007509437.tif 000233
【図３５】
2007509437.tif 000234
【図３６】
2007509437.tif 000235
【図３７】
2007509437.tif 000236
【図３８】
2007509437.tif 000237
【図３９】
2007509437.tif 000238
【図４０】
2007509437.tif 000239
【図４１】
2007509437.tif 000240
【図４２】
2007509437.tif 000241
【図４３】
2007509437.tif 000242
【図４４】
2007509437.tif 000243
【図４５】
2007509437.tif 000244
【図４６】
2007509437.tif 000245
【図４７】
2007509437.tif 000246
【図４８】
2007509437.tif 000247
【図４９】
2007509437.tif 000248
【図５０】
2007509437.tif 000249
【図５１】
2007509437.tif 000250
【図５２】
2007509437.tif 000251
【図５３】
2007509437.tif 000252
【図５４】
2007509437.tif 000253
【図５５】
2007509437.tif 000254
【図５６】
2007509437.tif 000255
【図５７】
2007509437.tif 000256
【図５８】
2007509437.tif 000257
【国際調査報告】
2007509437.tif 000258