JP2007509437A - 最適化システム - Google Patents

Abstract

【課題】一般的に適用できる効率的な最適化の手法を提供する。
【解決手段】本発明は、準最適解または最適解を求める問題を数学的にモデル化する高次元領域の関数を最小化するアプローチによって問題を最適化するシステムを提供する。このシステムは、コンピュータ・プログラムを実行するコンピュータ・システムと、最適化プログラムと、を有する。この最適化プログラムは、最適化問題を受け取り、連続型反復パラメータにより前記最適化問題をパラメータ化して連続最適化問題を生成し、連続最適化問題の解法を採用して前記最適化問題の最適解に向かう軌跡を定め、最適化問題の複数のパラメータを制御変数として定式化し、制御変数を計算するための反復法により最適な軌跡を含む所定領域内で準最適軌跡を求め、準最適軌跡を繰り返し計算して大域的な最適解を含む所定区間内の準最適解を得る。
【選択図】図２

Description

本発明は、さまざまな問題に対する相対的に最適な解を決定するため、および複雑な問題の解を求める計算プロセスと計算アプローチの構築を制御するために使用する一般最適化問題に関し、特に、一般的な目的関数のクラスを最小化するため、時間の経過につれて解を調整するため、およびさまざまな計算プロセスに対する相対的に最適な計算制御を構築し、実現するための一般的な計算の方法およびシステムに関する。

本出願は、いずれも２００２年１０月２３日に出願された仮出願第６０／４２０，９２１号および６０／４２０，９２２号の利益を主張する２００３年１０月２３日に出願された米国特許出願第１０／６９３，７２９号の一部継続出願である。

現在、数学および計算による最適化の分野には、膨大な研究および開発の労力が費やされている。最適化の技術は、複雑な系および組織向けの、費用対効果が高く、時間的に効率のよい制御戦略を含めて、さまざまな日常の問題および複雑な理論上の問題に適用され、こうした問題を解決している。たとえば、大規模なメーカー組織における労働力資源と設備資源のスケジューリングに最適化の技術を適用することによって、指定された量の製品を費用対効果が最も高い方法で生産し、さまざまな時点で計画された製品需要を満たすことができる。別の例として、都市の信号機制御のように大規模な制御問題に最適化の技術を適用することによって、都市全体にわたってさまざまな交差点での過度な遅延の発生をできるだけ抑えながら膨大な量の交通を処理できるようにし、さまざまな２次的および３次的な経路を利用して通行できるようにする。多くの自然体系は、物理的な特性と原理によって制御され、多くの物理的な制約によって制限された複雑な定義による問題の相対的な最適解を求めていると考えることができる。たとえば、複雑な化学反応は、分子などの個々の要素が最小のエネルギー状態と最大のエントロピーを追求するという比較的単純な熱力学の原理によって制御されており、分子が量子力学に従ってそれ自体を再構成し、それ自体の間でエネルギーを交換できるさまざまな方法によって制約されると考えることができる。さらに複雑な生体系では、化学と熱力学の原理だけでなく、情報量と情報交換、区分化とモジュール設計、および他に考慮すべき問題によっても制御されると考えることができる。このように、最適化問題には、自然現象や、複雑な系の設計、体系化、および制御や、情報の体系化、送信、および処理を、数学および計算によって表現した極めて多種多様なモデルが含まれる。

数理的にモデル化された問題の準最適解（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）によって導かれた軌跡を使用して得られた関数を繰り返しによって最適化し、超次元問題領域内にある１つ以上の局所的最小値を識別する。降下フィールド・メカニズムのさまざまな戦略を使用して、存在する可能性のある多数の局所的最小値を検出して評価し、数理的にモデル化された問題の準最適解を特定する。本発明のさまざまな実施形態では、さまざまな検索戦略を採用している。本発明の実施形態を表す方法でも、現在採用されている他の最適化法でも、指定された許容範囲内または妥当な時間内で、妥当な量のコンピュータ・リソースを使用して、広域的最小値を必ず検出することを保証しているわけではないことに留意するのは重要である。しかし、多くの実際的な最適化問題には境界が指定されており、局所的最小値の数は限られている。最終的に、指定された量のコンピュータ・リソースを使用して、指定された許容範囲または最大時間内で最善の解が見つかると、後処理のステップ２０４でこの解を変換し、ステップ２０１で初めに提供された問題に対する数理モデルのコンテクストで１つまたは複数の準最適解を提供する。

数理的なモデル化のステップと前処理のステップ、すなわち図２の２０１と２０２は、それぞれ具体的、仮説的な問題を使用しておそらく最もよく理解される。図３は、数理的なモデル化のステップと前処理のステップを説明するために使用する仮説的な問題を示している。この仮説的な問題には、ｐ個の製造工場３０２〜３０４、ｗ個の卸売店３０６〜３１１、ｃ個の顧客３１４〜３１８が存在する。ある製造工場、たとえば製造工場３０２は、図３のさまざまなタイプの線（たとえば製造工場３０２と卸売店３０６を結ぶ破線３２０）で示されるように、ｗ個の卸売店のいずれにも商品を出荷できる。さらに、各卸売店は、さまざまなタイプの線、たとえば卸売店３０６と顧客３１４を結ぶ破線３２２で示される卸売店と顧客の間のパスに従って、ｃ個の顧客のいずれにも商品を出荷できる。図３には、さまざまなスタイルの線を使用して、特定の製造工場からさまざまな卸売店へのパス、および特定の卸売店からさまざまな顧客へのパスが示されている。

図４〜７は、仮説的な問題のその他の側面と特性を、図３の描画法を使用して示している。たとえば、ｘ^ＰＷ _ｉｊの表記は製造工場ｉから卸売店ｊに出荷された商品の量を示し、ｘ^ＷＣ _ｊｋの形の項は卸売店ｊから顧客ｋに出荷された商品の量を示している。ｗ個の卸売店の候補が存在するが、ｗ個すべての卸売店を使用する必要はない。ブール変数ｂ_ｊは卸売店ｊが使用されているかいないかを示し、ブール値「１」すなわち「ＴＲＵＥ」は卸売店が使用されていることを、ブール値「０」すなわち「ＦＡＬＳＥ」は卸売店が使用されていないことを示す。図４では、たとえば、使用されている卸売店を示す４角形の実線の枠と使用されている卸売店に関連付けられられたブール値「１」で示されるように、卸売店１（３０６）、３（３０８）および５（３１０）が使用されており、卸売店２（３０７）、４（３０９）、および卸売店ｗ（３１１）を含めて卸売店５より後のすべての卸売店は使用されていない。図４では、使用されている卸売店は括弧内に連番で別の番号が付けられており、３つの使用されている卸売店には、順に１、２、およびｎの番号が付けられている。ただし、ｎは使用されている卸売店の数であり、図４の例では３である。

図５は、製造工場から卸売店に商品を出荷する際のコスト、および卸売店から顧客に商品を出荷する際のコストを示している。Ｃ^ＰＷ _ｉｊの形の項は製造工場ｉから卸売店ｊに商品を出荷する品目単位のコストを示し、Ｃ^ＷＣ _ｊｋの形の項は卸売店ｊから顧客ｋに商品を出荷する品目単位のコストを示している。図５では、出荷ルートのそれぞれについて、対応する品目単位のコストが表示されている。

図６は、仮説的な問題の付加的な特性を示している。各製造工場には製造能力が関連付けられており、製造工場ｉの出荷能力を示すｓ_ｉの形の項がある。各卸売店には受け入れ能力があり、これは卸売店ｊが商品の受け入れ能力を示す項ｖ_ｊで表される。また、各卸売店には運用コストがあり、卸売店ｊの運用コストを示すｒ_ｊの形の項がある。最後に、各顧客には商品の需要が関連付けられている。ｄ_ｋの形の項は、顧客ｋの商品に対する需要を示す。図６には、製造工場、卸売店、および顧客に関連する供給能力、受け入れ能力、運用コストおよび需要の詳細がすべて明確に示されている。最後に、図７は２つの製造工場、３つの潜在的な卸売店、５つの顧客からなる単純な仮説的問題を示しており、関連のコスト、能力および需要を数字で示している。

図３〜７に示す最適化問題は、顧客のすべての需要を満たし、能力が超えない範囲での卸売店の運用、製造工場から卸売店への出荷量、卸売店から顧客への出荷量に関する最小限のコスト構成を算出することである。最適化問題のより形式的な記号形式による数理的な表現は以下のとおりである。

まず、特定のシステム構成に関連する総コストは次の式で表される。

つまり、総コストは製造工場から卸売店、卸売店から顧客へのすべて出荷に関連するコストと、卸売店の運用コストの合計である。この仮説的な最適化問題において、総コストの最小値を求める。さらに、このモデルにはさまざまな制約が含まれる。まず、製造工場は製造工場の能力を超える量の商品を出荷することはできない。これは次のように表される。

すべての顧客は、少なくともその需要に等しい量の商品を受け取る。これは次のように表される。

さらに、製造工場から各卸売店に出荷する商品の量は、卸売店の商品受け入れ能力以下である。これは次のように表される。

同様に、卸売店から出荷する商品の量は、卸売店の商品受け入れ能力以下である。これは次のように表される。

卸売店には商品を保管しない。各卸売店は、受け取った商品をすべて出荷する。これは、次に示す保管に関する制約によって表される。

図３〜７に示す出荷のパスは単方向である。つまり、商品は返品されない。これは次のように表される。
ｘ^ＰＷ _ｉｊ≧０
ｘ^ＷＣ _ｊｋ≧０

仮説的な問題では、顧客の需要を満たすのに十分な製造能力と卸売店の受け入れ能力があることが仮定される。これは次のように表される。

商品は、運用している卸売店からのみ出荷できる。これは次の２つの制約によって表される。
ｘ^ＰＷ _ｉｊ(１−ｂ_ｊ)＝０
ｘ^ＷＣ _ｊｋ(１−ｂ_ｊ)＝０

運用している卸売店の数は、卸売店の候補の総数ｗ以下である。これは次のように表される。

最終的に、項ｙ_ｊはブール値（Boolean value）ｂ∈{０，１}である。

前述の数理モデルは、コストの関数と関連の制約をより簡潔に表すように若干変形してもよい。まず、出荷する商品の量を表すｘ^ＰＷ _ｉｊおよびｘ^ＷＣ _ｊｋのタイプのパラメータを、次のようにスケーリングされた出荷の項で置き換えてもよい。
ｍ^ＰＷ _ｉｊ＝ｘ^ＰＷ _ｉｊ／ｓ_ｉ、０≦ｍ^ＰＷ _ｉｊ≦１
ｍ^ＷＣ _ｊｋ＝ｘ^ＷＣ _ｊｋ／ｓ_ｋ、０≦ｍ^ＷＣ _ｊｋ≦１
前述の卸売店の能力に基づく制約は、冗長になるので制約のセットから削除でき、問題のモデルまたはモデルに関して得られる解に影響することはない。卸売店に商品を保管しないことを表す保管の制約は、スケーリングされた出荷の項を使用するために次のように変更される。

運用されていない卸売店からは商品を出荷しない。これは次の制約によって表される。
ｍ^ＷＣ _ｊｋ≦ｂ_ｊ
総コストは、商品のスケーリングされた出荷の項を使用して次のように表すことができる。

能力の制約は、商品のスケーリングされた出荷の項を使用して次のように表される。

前述の能力に関する制約は、次のように、０に関する不等式制約として作成し直すことができる。

ここで、運用している卸売店の数は、卸売店の候補の総数以下でなければならない、商品は返品できないこと、さらにｂ_ｊ項はブール値であるという事実は次のように表される。

図８は、図２に示す高レベルの制御流れ図のステップ２０２で呼び出される前処理のルーチンを示すより詳細な制御流れ図である。図８に示す前処理ルーチンの各ステップについては、前述の仮説的な数理モデルの対応する変換に関連付けながら以下で詳細に説明する。まず、ステップ８０１で、前処理ルーチンは、制約を含めて最適化問題を最初に数理的に表現した記号形式のモデルを受け取る。ステップ８０３で、前処理ルーチンは冗長な制約をできるだけ除去するように試みる。冗長な制約の除去の例は、上の変換されたモデルで実行されている。つまり、元の形のモデルにおいて冗長であると認識された卸売店の能力に関する制約は、変換されたモデルからは削除されている。次にステップ８０５で、すべてのブール値および整数値の決定変数（たとえば上の仮説的な例におけるブール変数ｂ_ｊ）は、浮動小数点変数に変換される。ステップ８０７で、結果として得られる数理的に表現された記号形式のモデルは標準の形式に配列される。ステップ８０９で、等式制約はゆるやかな不等式制約に置き換えられる。各等式制約は、正の値をとる追加の変数を使用する３つの不等式制約に緩和される。ステップ８１１で、追加の変数ｚと追加の制約が追加され、基準が追加の制約に変換される。最後にステップ８１３で、項を共に加えることによって、変換された最適化問題を内点法（Interior Point Method）のバリア関数（Barrier function）に変換する。各項は、制御変数をかけた制約を表す式に関するバリア関数の評価である。内点法のバリア関数には、２つの制御変数のセット、ｕとｒが含まれる。前処理ルーチンの結果として、数理モデルは、図２のステップ２０３に示す反復最適化の技術によって最適化できる形の記号形式で表現される。

図９は、図７に示す特定の最適化問題を最初に数理的に表現した記号形式のモデルを示している。コスト関数９０２は、整数変数ｘおよびブール変数ｂのベクトルに関して最適化される。制約には、製造工場の能力に関する制約９０４、顧客の需要に関する制約９０６、保管に関する制約９０８、最適化で選択されていない卸売店には出荷も入荷も発生しない、すなわち対応するブール変数はゼロであるという事実をモデル化した制約９１０、３つの卸売店のうち２つだけが運用されることを示す制約９１２、顧客から卸売店へも卸売店から製造工場へも商品の積み戻しがないことを示す制約９１４、およびｂ_ｊがブール値をとる変数であることを示す制約９１６が含まれる。

前処理ルーチンの第３のステップである図８のステップ８０５には、ブール変数、整数変数、およびその他の浮動小数点以外の変数を浮動小数点変数に変換するステップが含まれる。各ブール変数ｂ_ｉは、範囲［０，１］にあり、追加の制約を伴う浮動小数点の値ｙ_ｉに置き換えられる。
ｙ_ｉ（ｌ−ｙ_ｉ）＝０
追加の制約は、浮動小数点変数ｙ_ｉが０または１の２つの値のうち、いずれか１つしか取ることができないことを示している。各整数変数は、ｋ個の項からなる等価な２進多項式（base-2 polynomial）に置き換えられる。ここでｋは整数変数がとり得る最大値の２を底とする対数にほぼ等しい。２進多項式の各項は２のべき乗にブール値の係数をかけたもので構成される。したがって、２進多項式は［０，１］のすべての範囲で次に示す形のｋ個の制約を伴うｋ個の浮動小数点変数ｙ_０，ｙ_１，．．．ｙ_ｋ−１で置き換えられる。
ｙ_ｉ（ｌ−ｙ_ｉ）＝０
さまざまな型の決定変数の浮動小数点数への変換は、形式的に次のように要約される。

図１０は、仮説的な数理モデルから浮動小数点変数のみを使用したモデルへの変換を示している。

次に、数理モデルが標準の形式に配列される。標準の形式によるシンボリックな問題の表現は、形式的には次のように定義される。

つまり、標準の形式のモデルには、コストの関数Ｃ（ｘ）、Ｍ_１個の不等式制約

および複数個の等式制約

が含まれる。
図１１は標準の形式に変換した仮説的な問題の数理モデルを示している。元の等式制約と不等式制約を標準の形式で必要な右辺が０となる制約に変換するのは、まだ標準の形式になっていない制約式の両辺に項を加えるだけの簡単な作業であることに留意されたい。

前処理ルーチンの次のステップ（図８の８０９）は、すべての等式制約を不等式制約に変える操作である。これを行うために、各等式制約を２つの不等式制約に変換する。各不等式制約には、元の等式制約の左辺の式と正の束縛変数（binding variable）が含まれる。等式制約ｈ_ｊ（ｘ）から２つの不等式制約への変換は、形式的に次のように表される。

さらに、制約
ｓ_ｊ≧０
が追加され、ｈ_ｊ（ｘ）を表す制約が実現可能であることを保証する。

したがって、以下の等式制約から不等式制約への変換の後で、制約の総数は問題の標準の形式では元の等式制約の数に対して２ｋだけ増える。図１２Ａ〜Ｄは、変換された等式制約における束縛変数の効果を表している。図１２Ａは、値０を持つ元の等式制約ｈ_ｊ（ｘ）を示している（１２０２）。前述のように等式制約ｈ_ｊ（ｘ）が３つの不等式制約に変換されると、ｈ_ｊ（ｘ）の値の範囲は０から−ｓ_ｊ １２０４とｓ_ｊ １２０６の間にあるすべての値に拡張される。したがって、等式制約を３つの不等式制約に変換することによって、実質的に等式制約は点から区間に緩和される。後述するように、等式制約を緩和すると、等式制約によって１つまたは複数の局所的最小値に向かう最初の降下を過剰に拘束したり支配したりすることはなくなるので、反復最適化法が円滑化される。また、等式制約によって最適化の軌跡は問題領域の境界になることを余儀なくされるので、等式制約を緩和すると内点降下も促進される。しかし、最適化が収束し始めたら、等式制約の緩和を排除する必要がある。これは図１２Ｃおよび１２Ｄに示すように、バリア変数ｓ_ｊの絶対値を減少し、ｈ_ｊ（ｘ）の値に関する制約をさらに強めて０を中心にして対称的に配置される小さな領域に拘束することで達成できる。等式制約の緩和の程度は、基準に加えられた各束縛変数ｓ_ｊの係数ｒ_ｉで表すことができる。等式制約の不等式制約への変換は、形式的に次のように表される。

図１３は、仮説的な問題の等式制約を不等式制約に変換した後の数理的に表現された記号形式のモデルを示している。図１４は、束縛変数ｒ_ｉの効果を図形的に示している。図１４では、最適化法は初期点１４０２から始まり、カップ型の面１４０６で囲まれた問題領域１４０４内を最小１４０８に向かって降下する。最初は、束縛変数ｓ_ｉの係数ｒ_ｉは小さく、図１４において円盤形の領域１４１０で示されるように、初期点１４０２に関する値の範囲全体は大きい。問題領域を降下して最小１４０８に近づくにつれて、最適１４０８の近くでは、束縛変数ｓ_ｉの係数ｒ_ｉは大きくなり、解の集合に対する制約はしだいに厳格になる。最適化法が進むに従ってｓ_ｉはどんどん小さくなるので、変換された等式制約は再び極端に厳しく制約される（１４１２）。

前処理ルーチンの次のステップには、問題の数理モデルにエンベロープ変数（envelope variable）ｚと追加の制約を追加する処理が含まれる（図８の８１１）。エンベロープ変数ｚは、次のように定義される。

したがって、次の不等式制約が得られる。

エンベロープ変数ｚを導入した後、問題の形式的な記述は次のようになる。

図１５は、エンベロープ変数ｚを導入した後の仮説的な問題の数理モデルを示している。

図１６Ａ〜Ｄはバリア関数の概念を示している。たとえば、最適化問題の問題領域が図１６Ａに示すように放物線に基づくくさび形１６０２であると仮定する。最適化の初期点１６０４では相対的に大きなｚ値１６０６をとり、最小１６０８で０に減少する。図１６Ｂは、問題領域１６０２のｘ−ｙ平面への投影図１６１０を示している。初期点１６０４から最小点１６０８までの最適化の軌跡は３つの線分１６１２〜１６１４で示されている。内点法では、問題領域を移動するとき、最適化の軌跡は問題領域内にとどまる必要がある。換言すると、問題領域内で最小に至る降下の軌跡上にある現在の点は、常に問題領域内すなわち超次元問題領域内の多様体の境界の内部に残っている必要がある。たとえば、図１６Ｂで示されるように、最初の軌跡１６１２は問題領域に対応する形状の境界に向かう角度で進み、点１６１６で境界と交わる。引き続き同じ方向に進むと、境界を越えて破線の線分１６１８に沿って問題領域外のある点を通過することになる。最適化の軌跡を問題領域内に維持するには、問題領域内では影響が小さいか全くないが、最適化の軌跡が問題領域の境界面の近くになると、軌跡を境界から問題領域に戻すために値が増大するようなバリア関数を問題のモデルに含める必要がある。たとえば図１６Ｂでは、最適化の軌跡１６１２が点１６１６に近づくに従って、数理モデルに追加されたバリア関数が急速に増大し始め、エンベロープ変数ｚの値を増大することによって最適化の軌跡を問題領域の内側の境界面から離れるように屈曲させる。屈曲した後は、後続の最適化によって境界点１６２０に向かう新しい第２の軌跡に続く最適化の軌跡１６１３が確立される。境界点１６２０の直前でバリア制約によって再び目的関数ｚの値が大きくなって再び最適化の軌跡の向きが問題領域の内側に屈曲し、新しい最後の軌跡１６１４を得る。

図１６Ｃ〜Ｄはバリア関数の動作を示している。図１６Ｃに示すように現在の最適化の軌跡１６２２を想定し、図１６Ｄに示すように、図１６Ｃに示す現在の軌跡１６２２に沿って問題の領域内を移動する目的関数ｚのｘ軸に関する値について考察する。図１６Ｄからわかるように、ｚの値は現在の軌跡１６２２に沿った問題領域の内部では比較的小さい。しかし、ｘが減少して左の境界上の点１６２４に近づくか、増大して右の境界上の点１６２６に近づくと、バリア関数が影響してそれぞれのｚの値１６２８と１６３０が急速に増大する。換言すると、軌跡を問題領域の内部に拘束するために、モデルに含まれる不等式のそれぞれについて、内点の目的関数にバリア関数を含む項が追加される。軌跡が問題領域から出た場合、つまり問題領域を囲む面の外側に達した場合は、本発明のいくつかの実施形態を表す最適化法は失敗する。

モデルにバリア関数に含む項を追加するのは、最適化の軌跡が問題領域内に確実にとどまるようにするための簡単な方法である。問題領域を囲む面を指定するのは制約なので、制約ごとにバリア関数が追加される。本発明のさまざまな実施形態で使用する最適化技術に適したバリア関数項は、こうした項に関連付けられた制御変数である正の変数の係数ｕ_ｉに制約の式ｇ_ｉの対数をかけ、負号を付けたものである。つまり、次のように表される。
−ｕ_ｉlog(ｇ_ｉ)
制約の式ｇ_ｉの値が０に向かって減少するに従って、制約の対数log(ｇ_ｉ)の値は急速に減少し、これに対応して絶対値は急速に増大する。Log(ｇ_ｉ)の負号を付けることで、制約式ｇ_ｉが０に非常に近くなると、急速に増大する値が得られ、極限では無限大に近づく。したがって、制約がさらに減少して０になるのを防止できる。この変数の係数ｕ_ｉは、前述の束縛変数の係数ｒ_ｉと同様に、制御変数とみなすことができる。最適化が局所的最小値に向かって最も効率のよい軌跡をたどるように、最適化の間に係数ｕ_ｉとｒ_ｉが変更される。しかし、係数ｕ_ｉは係数ｒ_ｉと違って初期値が大きく、本発明の１つの実施形態によって最適化の軌跡が局所的最小値に近づくに従って減少する。これは、局所的最小値はしばしば問題領域の境界面に近いか境界面上にあるためである。目的関数にバリア関数を追加する形式的な式は次のようになる。

図１７は、バリアの制約を追加した後の仮説的な問題の記号形式のモデルを示している。バリア関数を追加することによって、最初に説明した実施形態の前処理ルーチンを具体化した前処理のステップは完了する。

次に、本発明の１つの実施形態を表す中心的な最適化法について概説する。図１８は、問題領域内の広域的勾配降下フィールドを使用して、初期点から局所的最小値に至る最適化の軌跡の方向を示すようすを表している。問題領域内の各点、たとえば点１８０２などで、目的関数Ｆ(ｘ，ｓ，ｚ，ｕ，ｒ)の偏導関数ベクトルが、目的関数の決定変数ｘ_ｉ、束縛変数ｓ_ｉおよびエンベロープ変数ｚに関する１次偏導関数のベクトルとして計算される。同様に、決定変数ｘ_ｉ、束縛変数ｓ_ｉ、およびエンベロープ変数ｚに関するＦ（ｘ，ｓ，ｚ，ｕ，ｒ）の２次偏導関数の行列として与えられる、関数Ｆ（ｘ，ｓ，ｚ，ｕ，ｒ）のヘッシアン（Hessian）行列は、領域のすべての点で計算可能である。目的関数の勾配は、領域内の各点においてヘッシアン行列の逆行列とＦ（ｘ，ｓ，ｚ，ｕ，ｒ）の偏導関数の積で得られるベクトルとして定義される。これは、決定変数と束縛変数に関して目的関数の値が降下する方向を表しており、点１８０２からの方向１８０４を示している。勾配は問題領域内の各点で計算できるので、目的関数の勾配は問題領域内における降下方向を示し、最小へのパスを指定するベクトル場を表している。一般に、図１８の点１８０２のように局所的最小値から離れた点の勾配は絶対値が相対的に大きく、点１８０６のように局所的最小値に近い点の勾配は絶対値が相対的に小さくなる。

最適化法の数理的記述は、多くの表記法によって簡便化される。まず、目的関数Ｆ（ｘ，ｓ，ｚ，ｕ，ｒ）は、Ｆ（ｙ，ｕ，ｒ）として表現できる。ここで、

である。

したがって、目的関数は、決定変数、束縛変数、エンベロープ変数からなるベクトルｙ、バリア関数の係数からなるベクトルｕ、および束縛変数の係数からなるベクトルｒの関数と考えることができる。ベクトルｙは「状態ベクトル（state vector）」と呼ばれ、ベクトルｕとｒは「制御ベクトル（control vectors）」と呼ばれる。ベクトルｙはＮ個の要素を持っている。ただし、Ｎ＝ｎ（決定変数の数）＋ｋ（束縛変数の数）＋１（エンベロープ変数の数）である。

図１９は、最適化の軌跡を示している。軌跡１９０２は、ベクトルｒとｕの長さを横軸、目的関数Ｆの値を縦軸としてプロットされている。最適化の軌跡は、初期点１９０４から開始され、局所的最小値に向かって目的関数Ｆの値が降下する。降下する間に、ベクトルｕの長さで表されるバリア関数の値は０に向かって減少し、ベクトルｒの長さで表される束縛変数の係数の値は増加する。軌跡は長さΔτの線分に分割される（たとえば最初の線分１９０６）。最適化の間に、ベクトルｙとｕは目的関数の値を小さくするように継続的に再計算される。図１９では、ベクトルｒで表される束縛変数の係数の値も、各区間Δτの間に再計算されることが示されている。しかし、ここで説明する本発明の実施形態では、束縛変数の係数は最適化の軌跡の終わり近くでのみ調整される。

離散的な最適化の技術は次のように表すことができる。

ｙ_ｋ＋ｌ＝ｙ_ｋ＋φ_ｋ（ｙ_ｋ，ｒ_ｋ，ｕ_ｋ）

ここで、ｋ＝０，１，．．．，ｋ_convergenaeである。換言すると、計算による最適化法の繰り返しごとに、現在の状態ベクトルの値ｙ_ｋと、ｙ_ｋ、ｒ_ｋ、およびｕ_ｋに依存する離散関数φ_ｋとから新しい状態ベクトルの値ｙ_ｋ＋１が計算される。ただし、反復回数ｋは０からｋ_convergenceまでの値をとり、ｋ_convergenceで目的関数が局所的最小値となる点の状態ベクトルを決定できる。最適化の最初のステップでは、状態ベクトルと制御ベクトルｙ、ｕ、およびｒは連続型反復パラメータ（continuous iteration parameter）τを媒介変数として使用する。ここでτの範囲は０から反復限界（iteration horizon）Ｔまでである。連続型反復変数τをベクトルｙ、ｕ、およびｒの媒介変数として使用することにより、前述の繰り返しはｙ変数が従属変数であり、ｕとｒは制御変数である１次元微分方程式に変換される。微分方程式に基づく数理的な周知の技術を使用すると、この微分方程式の解として効率的な最適化の軌跡を計算できる。このことについては、以下で説明する。これで、目的関数は形式的に次のように表される。
Ｆ(ｙ(τ)，ｕ(τ)，ｒ(τ))
ここで、０≦τ≦Ｔである。たとえば、ｙ_ｋ＋１＝ｙ_ｋ＋φ_ｋ(ｙ_ｋ，ｒ_ｋ，ｕ_ｋ)がニュートン法による反復（Newton iteration）の場合は、対応する微分方程式は次の形である。

ニュートン法（Newton’s mothod）は、ｄｙ／ｄτの式を提供する微分方程式を構築するための方法の１つにすぎないことに留意されたい。この方法は、計算を単純にし、容易にするために選択されている。ただし、本発明の実施形態では、ｄｙ／ｄτを近似するために、他の多くのアプローチのいずれを使用してもよい。

以下の説明では、目的関数Ｆの偏導関数ベクトルＧとヘッシアン行列Ｈを頻繁に使用する。偏導関数ベクトルＧは、次のように定義される。

ヘッシアン行列ＨはサイズＮｘＮの行列であり、次のように定義される。

連続する反復変数τに関する状態ベクトルｙの正規化された変化の率は、ヘッシアン行列ＨとベクトルＧを使用して次のように表現される。

または、次のようにも表される。

ここで、ｗ(τ)＝Ｈ^−１Ｇ(ｙ(τ), ｕ(τ), ｒ(τ))、および
ρ(ｙ(τ), ｕ(τ), ｒ(τ))＝Ｇ^ＴＨ^−１Ｇ(ｙ(τ), ｕ(τ), ｒ(τ))
である。

最適化問題は、一般的な最適制御問題の特殊な事例である最短時間問題に書き換えることができる。換言すると、目的関数Ｆ（ｙ，ｕ，ｒ）を連続にすることによって、連続型方法を使用して、目的関数Ｆ（ｙ，ｕ，ｒ）の値を初期点から最小に導く状態変数の変化で構成される最適化の軌跡ｙ（τ）を計算するための計算のアプローチが得られる。しかし、問題領域が非常に大きく、超次元であるために、詳細な計算上の制約またはメタレベルの制約がないと、問題領域内の降下勾配フィールドの方法が得られても、計算による勾配フィールドの降下には収束が遅いという問題がある。

最適制御問題を定式化する。このことによって、収束限界を最小化し、制御変数に関する要件を満足することによって、収束を速めようとする。目的関数Ｆ（ｙ（τ），ｕ（τ），ｒ（τ））の最小化を制御するために、最適制御の定式化は、制御変数ｕ（τ），ｒ（τ）を使用し、収束限界を最小化する制御変数の値を選択することによって、最適化の軌跡を求める。最適化の軌跡を制御する最適制御問題は、ｕ（τ），ｒ（τ），τ∈［０，Ｔ］に関して形式的に次のように表現される。

ここで、Ｑは対称正定値行列（symmetric positive definite matrix）であり、最適制御問題は、次の制約条件に従う。

最小収束限界問題の目的関数Ｊ(ｕ(τ)，ｒ(τ)，Ｔ)の最小化により、Ｆ(ｙ(τ)，ｕ(τ)，ｒ(τ))を強制的に最小に収束する。この場合、最小の計算限界ＴにおいてＧ＝０であり、収束区間（convergence interval）は連続する反復変数の積分である。最適制御問題の解は、状態ベクトル変化（state-vector-change）の最適化の軌跡ｙ^＊(τ)を生成するための最適制御ポリシーｕ^＊(τ)およびｒ^＊(τ)、τ∈[０，Ｔ]を提供する。最小収束限界問題の目的関数Ｊ(ｕ(τ)，ｒ(τ)，Ｔ)を定義する積分の最初の項によって、最適化の軌跡は連続型反復変数τに関してできるだけ効率的な軌跡ｙ^＊に強制的に導かれる。以下の説明では、目的関数Ｆの偏導関数ベクトルＧとヘッシアン行列Ｈを頻繁に使用する。積分の第２項によって、ベクトルＧは強制的に０に向けられ、そのことによって目的関数Ｆは強制的に最小に導かれる。したがって、これらの項の組合せによって最適制御ポリシーｕ^＊(τ)とｒ^＊(τ)が生成され、このポリシーによって目的関数Ｆは最適な状態ベクトル変化の軌跡ｙ^＊（τ）をたどって強制的に局所的最小値に導かれる。ここで、最適制御問題（１）〜（５）について説明する。第１の制約は正規化された状態ベクトル変化の最も急な降下を表しており、元の最適化問題が新しい最短時間制御問題に組み込まれる基本的な手段である。第２と第３の制約は、制御ベクトルの値ｕが常に非負であり、特定の極大制御値ｕ_ｍａｘｉ未満であることと、ｕベクトルの値が最適化処理中に増大しないことを保証する。第４と第５の制約は、第２の内部的な制約と同様であり、束縛変数の係数ｒ_ｉに関連する。第４と第５の制約は、束縛変数の係数ｒ_ｉが最小値ｒ_ｍｉｎｉを下回らないことと、束縛変数の係数ｒ_ｉが最適化処理中に減少しないことを保証している。

最適制御理論の主な結果は、ポントリャーギンの最小原理（The Pontryagin Minimum Principle）として知られており、最小収束限界問題のフレームワークを提供し、状態ベクトル変化の最適化の軌跡ｙ^＊（τ）を検出するための最適制御ポリシーｕ^＊（τ）およびｒ^＊（τ）、τ∈［０，Ｔ］を生成する。このアプローチには、ハミルトニアン（Hamiltonian）関数が含まれる。この関数は次のように定義される。

ハミルトニアンは、ｕ^＊（τ）およびｒ^＊（τ）におけるｕおよびｒの関数としての最小値をとる。換言すると、ハミルトニアンは最適制御ポリシーの最小値をとる。制御ポリシーｕ^＊（τ）が最適な場合は、共役（costate）と呼ばれる絶対的に連続な関数Ｐ（τ）が存在する。このとき、状態関数ｙ(τ)およびと共役関数ｐ(τ)の導関数は、τ∈[０，Ｔ]となるすべてのτについて以下のように表され、

以下の境界条件

を満たしている。制御問題の解は、次のように要約される。

このハミルトニアンは、機械系の総エネルギーを表現する機械系のハミルトニアンに類似している。関数ｙ^＊（τ）は、連続型時間区間において機械系のコンポーネントのエネルギーが最小になる位置の指定に類似しており、共役関数ｐ^＊（τ）は、連続型時間区間において機械系最小のコンポーネントのエネルギー・モーメントの指定に類似している。上で導かれた結果を使用すると、最適化の軌跡に沿った偏導関数∂ｙ^＊／∂τの順方向（forward direction）の計算が可能である。偏導関数∂ｐ／∂τは、順方向の計算結果ｙ^＊（τ）から、τ＝Ｔによって後退型で計算できる。もう一度伝統的な機械の問題に例えると、機械系のコンポーネントのモーメントは、本質的にコンポーネントの時間的な位置の変化によって決定する。したがって、モーメントまたは共役は、順方向で計算した位置から後退型で計算する必要がある。

ハミルトニアンの最小化は、それ自体が最適化問題であり、ハミルトニアンが最小化される場合の最適な制御ポリシーを生成する。最小化問題は、ニュートン近似法を使用して次のように解くこともできる。

したがって、最短時間問題および最適制御問題を定式化することにより、連続する反復変数τの０から限界Ｔまでの範囲にわたって、状態ベクトルｙおよび制御ベクトルｕとｒを順方向で計算し、共役ベクトルｐを後退型で計算する方法が得られる。さらに、この操作は１つの反復最適化法となり、ここで目的関数Ｆは、最適化制御ポリシーｕ^＊（τ）およびｒ^＊（τ）の下で最適な状態ベクトルの軌跡ｙ^＊（τ）を経由して強制的に局所的最小値に導かれる。

図２０は、上記で導かれた状態ベクトル変化の最適化の軌跡ｙ^＊および最適化制御ポリシーｕ^＊とｒ^＊に基づいた反復最適化法を図式的に説明している。反復最適化法は、２つのネストしたループ、すなわち反復変数Ｋで制御される外側ループと反復変数ｋで制御される内側ループからなると考えることができる。外側のネストしたループは、目的関数Ｆを最小化する状態ベクトルｙを検出するための、連続する繰り返しの範囲０からＴまでの状態ベクトル軌跡および制御ポリシーの偏微分の積分を表している。このような外側ループの繰り返しごとに、連続する反復変数τの短い区間に関して制御ポリシーベクトルｕとｒおよび状態ベクトルｙに関する複数の順方向の計算がステップτ_ｋで実行される。ただし、ｋは０から最後のｋの値であり、換言すると、τ_０、τ_１、τ_２、．．．、τ_{ｆｉｎａｌ}をとり、τ_{ｆｉｎａｌ}≧Ｔになると、ｙ（τ_ｋ）とｕ（τ_ｋ）の計算結果を使用して、τ＝Ｔからτ＝０までの共役ベクトルｐ（τ_ｋ）を後退型で計算する。内側ループを実行するごとに、ｙ（τ_ｋ）の値の計算は前のｕ（τ_ｋ）とｒ（τ_ｋ）の計算に基づいて行われ、ｙ（τ_ｋ）とｕ（τ_ｋ）の両方の計算は直前の外側ループの繰り返しが終了したときに計算された共役ベクトルｐ（τ_ｋ）を使用して行われる。最初の繰り返しでは、制御共役値は定数で近似される。

図２０において、外側ループは外側の円２００２で表され、内側ループの計算は、たとえばステップ２００４などの階段状の放射線と円周の断片で表される。繰り返しの最適化プロセスは、状態ベクトル、制御ベクトル、共役ベクトルの初期値ｙ_０、ｕ_０、ｒ_０およびｐ_０（それぞれ２００６）で開始される。外側ループ変数Ｋは０に設定される（２００８）。次に、複数の階段状の内側ループの繰り返しがｋ＝０で開始され、内側ループを次々に繰り返すたびに内側ループ変数ｋが１ずつ増分される。最初の内側ループの繰り返しで、制御ベクトルｕ（τ_ｋ＋１）およびｒ（τ_ｋ＋１）が計算される（２００４）。さらに、τ＝τ_ｋ＋１、ｙ（τ_ｋ＋１）における状態ベクトルの値が多くの連続型計算２０１０〜２０１６で計算される。最初の内側ループの繰り返しが終わるときに、局所的最小値に到達したかどうかが決定される。達した場合は反復最適化法が終了する。この様子は、図２０の破線２０１８で示されている。達していない場合は、次の内側ループの繰り返し２０２０〜２０２１が行われる。内側ループの繰り返しは、τ_ｋ≧Ｔ（２０２２）となるまで続く。この時点で共役ベクトル値ｐ（τ）が、τ＝Ｔ からτ＝０までの複数のステップ２０２４−２０３０で計算され、外側ループの繰り返し制御変数Ｋが増分されて完全な外側ループの繰り返しが終わり、次の外側ループの繰り返しを開始する点２０３２に戻る。

図１９と２０を併せて考察することによって、繰り返しによる最適化プロセスを図式的、概念的に見ることができる。最適化プロセスの内部ループの繰り返しは、そのそれぞれが図１９に示されている最適化の軌跡における一部分、たとえば一部分１９０６に対応すると考えることができる。外側ループ全体の１回の繰り返し（τ＝０からτ＝Ｔまで）によって、最適化の軌跡の大きな一部分が生成されると考えることができる（たとえば初期点１９０４から始まって点１９０８まで下方に伸びる最適化の軌跡の一部分）。したがって、複数の外側ループの繰り返しによって、複数の大きな最適化の軌跡の一部分を計算する。この一部分は、最終的に状態ベクトルｙを目的関数Ｆが少なくとも局所的最小値となる値に向かわせる。図２１は、反復最適化法全体を示す追加の図式的な説明を提供する。図２１に示されているように、初回の外側ループの繰り返しは、初期点２１０２から点２１０４までの最適化の軌跡の一部分を生成する。２回目の外側ループの繰り返しは、点２１０４から点２１０８までの最適化の軌跡の一部分２１０６を計算する。３回目の外側ループの繰り返しは、点２１０８から点２１１２に至る最適化の軌跡の一部分２１１０を生成する。限界の値Ｔは各外側ループの繰り返しごとに変化してもよいことに留意されたい。このように、図２１では、各最適化の軌跡一部分の終点はＴ_ｋ（ただし、Ｋ＝１、２、および３）とマークされており、外側ループによる繰り返しの終点はそれぞれ次回の外側ループによる繰り返しの開始点となることが示されている。各内側ループのステップにおけるΔτのサイズも変化してよいことにも留意されたい。

図２２は、一般的な反復最適化のステップ（図２の２０３）を説明する制御流れ図である。まず、ｙ、ｕ、ｒおよびｐの初期値が選択され、外側ループの反復変数Ｋが初期値０に設定され、さらに初期限界Ｔ_０がステップ２２０２で計算される。次に、内側ループを形成するステップ２２０４〜２２０８で、τ＝０からτ_ｋずつ増分されてτ＝Ｔ_ｋに至る間に、制御ベクトルｕとｒおよび状態ベクトルｙは順方向で伝達される。まず、ステップ２２０４で制御ベクトルｕとｒの次の値が計算される。さらに、ステップ２２０５で、τ＝τ_ｋからτ＝τ_ｋ＋Δτ_ｋに状態ベクトルｙが伝達される。ステップ２２０６で最適化が収束したがどうかが判定される。収束した場合は、反復最適化が完了したので復帰する。収束していない場合は、ステップ２２０７でτ_ｋ＋Δτ_ｋがＴ_Ｋより小さいかどうかが判定される。そうでない場合は、ステップ２２０８でτ_ｋ＋ｌにτ_ｋ＋Δτ_ｋが代入され、内側ループの制御変数ｋが増分される。しかし、τ_ｋがＴ_Ｋ以上の場合は、ステップ２２０９で共役の値がＴ_Ｋから０まで逆順に計算され、ステップ２２１０で、反復変数および状態ベクトル、制御ベクトルの計算された値を含むデータベースが更新される。

凸問題（convex problems）では臨界点は１つしかないので、図１９〜２１を参照しながら以上で説明した反復最適化法は、数千もの決定変数を持つ極めて大きな超次元問題領域に対しても、通常は妥当な時間内に唯一の軌道で収束する。しかし、適用例における最適化問題は通常凸ではなく、問題領域は多くの臨界点を持つという特徴がある。図２３は、このような問題領域を示している。図２３において、問題領域２３０３の境界面は、複数の中央が突出した突起部（たとえば突起部２３０４）を持つ椀状の形である。これらの突起部は、３つの頂点（２３０４、２３０６、２３０８）がある小さな山岳地帯と、深い局所的な谷、すなわちくぼみ２３１０と考えることができる。問題領域の大域的な最適化の価値は、外側のくぼみの真中にある山岳地帯の下方深くに存在する。

図２４は、図２３と同じ問題領域を示しているが、臨界点が明示的に示されている。山岳地帯２３０４、２３０６および２３０８の頂上は問題領域内の極大を表しており、小さな谷の底２３１０は局所的最小値を表している。山岳地帯内の鞍点（たとえば鞍点２３１２および２３１４）も臨界点である。図２５は、臨界点（たとえば図２４における臨界点２３０４）の拡大表示である。図２５に示すように、臨界点２５０６に近いいくつかの点２５０２〜２５０４で計算されたｘ方向およびｙ方向の偏微分はすべて絶対値が小さい。したがって、臨界点２５０６に近い点２２０２〜２３０４で計算された勾配も０に近い。要約すると、τに関する目的関数の導関数

によって、臨界点に近い最適化の軌跡を決定する手段が得られる。したがって、点が臨界点に近づくにつれて勾配が０になるかまたは０に近づくと、τに関する状態ベクトルの軌跡の変化（ｒａｔｅ）も０となる。さらに、Ｇは臨界点の近くでは敏感であり、したがって小さな局所的な変化がＧの比較的大きな変化の原因となることがあり、Ｇの符号が変化する場合もある。しかし、Ｇが小さい場合は、状態ベクトルの軌跡を計算するために使用するニュートン法や他の同様の方法は、一般的に最適化の軌跡を発見するためには有効ではない。換言すると、勾配フィールドは一般的に準最適軌跡または最適な軌跡の方向を提供する。したがって、勾配フィールドが０の場合は、可能なすべての軌跡が同等に適切であるように思われる。これは、磁北極または磁南極でコンパスを使おうとする場合と同様の状況である。後で説明するように、臨界点の近くでは別の方法を使用して状態ベクトルの最適化の軌跡を計算する必要がある。

臨界点については、別の考慮すべき問題がある。この考慮すべき問題には、問題領域におけるいくつかの可能性のある局所的最小値の識別と評価を試みるために、大域的な最適化法で臨界点の近くで使用する必要がある検索に似た戦略が含まれる。図２６は、第１の検索技術を示している。図２６では、水平線Ｔｎへ向かう本最適化法の外側ループの繰り返しは状態ベクトルを点２６０２に伝達する。この点は、山岳地帯の頂上である臨界点２３０６の近くである。この点における勾配が０に近いという事実によって、最適化法に対して臨界点の近くに到達したことを伝える。ただし、次のステップは到着した臨界点のタイプによって決まる。たとえば点２３０６のように極大の場合は、本最適化法では何らかの方法を使用して最適化の軌跡を大きく変更し、極大２３０６を離れて局所的最小値に向かう降下２３０８を続ける必要がある。しかし、たとえば局所（local）２３１０のような局所的最小値の近傍では、本最適化法は局所的最小値２３１２から別の点２３１４に上昇し、この点から局所的最小値に向かって降下２３１６を続けるような技術を使用する必要がある。

図２７は、問題領域を検索するための第２の戦略を示している。図２７では、検索方法は極大２３０６に近い点２６０２に再度到着している。しかし、極大を検出すると、本最適化法は異なる方向２７０２〜２７０５に向かって降下を続ける相対的に独立した最適化の軌道をいくつか追加して生成する。相対的に独立した軌道の１つである２７０４は、鞍点２３１４に近い点２７０８に到達する。この点で、この軌道２７０４は追加の軌道を分岐または生成し、結果として得られる相対的に独立した各軌道は鞍点から異なる方向２７１０と２７１２に向かって降下する。最初に生成された別の相対的に独立した軌道２７０２は、局所的最小値２３１０に到達する。この時点で、新しい最適化の軌道を生成せず、相対的に独立した軌道２７０２は終了し、親軌道２７０１に局所的最小値を報告する。親軌道は子によって識別された最小の局所的最小値を選択し、それを局所的最小値としてその親に返す。結果として問題領域の再帰的な検索が行われ、最適化問題の準最適解または最適解が識別される。別の実施形態では、親軌道は報告された局所的最小値の情報を使用して問題領域を切り詰め、局所的最小値に基づいて問題領域を修正でき、問題領域が切り詰められたことを相対的に独立した他のすべての軌道に伝えることができる。

本最適化法は、図２６に示す唯一の軌道の技術でも、図２７に示す複数の軌道の技術でも、勾配が０に近い値に減少したときに到達した臨界点の種類を決定する手段が必要である。臨界点の種類を解明するための臨界点指標（critical-point beacon）は前述のヘッシアン行列である。勾配が０に近づいたときに、ヘッシアン行列のすべての要素が負ならば、極大に到達している。勾配が０で、ヘッシアン行列のすべての要素が正ならば、局所的最小値に到達している。それ以外の場合に、勾配が０でヘッシアン行列の要素の符号が混在していれば、変曲点または鞍点に到達している。正則（non-singula）ヘッシアン行列の臨界点指標属性は形式的に次のように記述される。

点ｙについて、

のとき、
Ｈの全ての固有値が負ならば、ｙは局所的または大域的な最大値の近傍にあり、
Ｈの全ての固有値が正ならば、ｙは局所的または大域的な最小値の近傍にあり、
その他の場合、ｙは鞍点（saddle point）にある。

図２６と２７に関連して前述したように、唯一の最適化の軌道は、階層的に体系化されたモジュールの集合で構成される。このようなモジュールを離散型の計算環境で同時に実行し、前述の最適化法による状態ベクトルの伝達に関連する大量の計算タスクを分散してもよい。図２８は、１つの最適化の軌道を構成するモジュールとモジュール化された編成を示している。本明細書では、最上位レベルのモジュールを「スーパー（super）」２８０２と呼ぶ。スーパーには、中間結果の比較的大きなデータベース２８０４が含まれる。スーパーは、それぞれの連続する限界Ｔ_Ｋの計算、図２７を参照しながらすでに説明したように、初期点の選択（複数の軌道による最適化法におけるルート軌道の場合）、および下位モジュールの起動とアクティビティ調整を担当する。第１レベルの下位モジュールには、ｕ＿ａｎｄ＿ｒモジュール２８０６、ｙモジュール２８０８およびｐモジュール２８１０が含まれる。ｕ＿ａｎｄ＿ｒモジュール２８０６は、図２２のステップ２２０４に示すように、制御ベクトルｕとｒの順方向の計算を担当する。

ｙモジュール２８０８は、図２２のステップ２２０５に示すように、連続する区間Δτを通じた状態ベクトルの伝達を担当する。ｐモジュール２８１０は、図２２のステップ２２０９に示すように、反復最適化法の外側ループの繰り返しが終了するたびに、Ｔ_ｋから０までの後退型による共役ベクトル値の計算を担当する。

ｕ＿ａｎｄ＿ｒモジュール２８０６は、制御ベクトルの順方向の計算を、ｍ個のエージェント２８１２〜２８１５に分散している。ｙモジュール２８０８とｐモジュール２８１０は各モジュールのタスクをそれぞれｎ個のエージェント２８１６〜２８１８および２８２０〜２８２３に分散している。ｙモジュール２８０８がｙの伝達をｎ個のｙモジュール・エージェント２８１６〜２８１８に分散することは、最適化計算全体を分散し、パラレルな計算の機会を提供するために重要なアプローチである。ｙモジュール２８０８およびそのエージェント２８１６〜２８１８で行われる状態ベクトル伝達の計算は重要であり、最も計算コストがかかり、状態ベクトル変化による最適化の軌跡の計算に必要なタスクである。本発明の別の実施形態では、ｐモジュールは効率的な補間法を使用して共役ベクトルの値を推定するので、問題をｎ個のｐモジュール・エージェントに分散する必要はない。したがって、ｐモジュールは分散せずに共役を計算する。さらに別の実施形態では、ｐモジュールは共役の補間を少数のｐモジュール・エージェントに分散できる。

スーパー２８０２で保守するデータベース２８０４は大きいため、メモリだけでなく、軌道のモジュールを分散する１つまたは複数のコンピュータ・システムに直接接続またはリモート接続されたマス・ストレージ・デバイスに格納してもよい。図２９は、本発明の様々な実施形態においてスーパーが保守するデータベースの内容の一部を示している。一般に、スーパーは、０からＴまでの水平線の各部分区間τ_ｋで、モジュールによって計算されるさまざまな量の値を保守する必要がある。したがって、図２９では、たとえばベクトルＧのデータ構造２９０４の部分列（sub-column）２９０２など、データ構造の第１列の各部分列は、特定のτの値に対応するエンティティの計算値を表している。図２９において、部分列２９０２はτ＝０におけるベクトルＧの値を表している。スーパーは、状態ベクトルの値２９０６、制御ベクトルｕとｒの値２９０８〜２９０９、ベクトルｗの値２９１０、スカラーρの値２９１１、モース（Morse）フラグ・ベクトルの値２９１２、共役ベクトルｐの値２９１４、目的関数Ｆのスカラー値２９１５、ベクトルＧの値２９０４、ハミルトニアンの値２９１６、現在のＴと反復制御変数および他のそのような局所的な制御変数の値２９１８、およびヘッシアン行列の逆行列の値２９２０〜２９２７を保存する。

図３０〜３４は、図２８の描画法を使用して、本最適化法の外側ループが繰り返されるときに、階層的に構造化された最適化の軌道のモジュール間で情報が流れる様子を示している。ルート軌道の場合は、図３０に示すように、最適化の軌跡の初期点が決定する（３００２）。さらに、他のデータ構造の中で、状態ベクトル、共役ベクトル、制御ベクトル、および制御変数の初期値が決定され、さまざまなデータ構造と要素３００４が初期化される。ルート以外の軌道は初期点および、状態ベクトルや制御ベクトルなどのさまざまな計算エンティティの現在の値を親軌道から継承する。図３１に示すように、内側ループの最初のステップ（図２２の２２０４）で、ｕ＿ａｎｄ＿ｒモジュール２８０６はスーパー ３１０２から状態ベクトル、共役ベクトル、ベクトルｗおよびヘッシアンの逆行列の値を受け取り、これらの値および／またはこれらの値の一部をｕ＿ａｎｄ＿ｒモジュール・エージェント３１０４〜３１０７に分散し、中間の計算結果をｕ＿ａｎｄ＿ｒモジュール・エージェント３１０８、３１０９、３１１０、および３１１１から受け取り、新しい制御ベクトルの値を計算し、さらに新しい制御ベクトルの値３１１２をスーパーに返す。図３２は、内側ループの次のステップ（図２２の２２０５）におけるデータの流れを説明している。このステップでは、ｙモジュール２８０８はスーパー ３２０２から新しく計算された制御ベクトルｕとｒを受け取り、これらの値にヘッシアンや他のベクトルなどの値を含むスーパーから取得したその他の値および／またはｙモジュールが計算したその他の値を追加してｙエージェントに配信する（３２０４〜３２０６）。次に、中間の計算値をｙエージェントから受け取り（３２０８、３２０９、３２１０）、特に、状態ベクトル、ｗベクトル、ｗベクトルの状態ベクトルに関する微分、ヘッシアンの逆行列を計算してスーパーに返す（３２１０）。ｙモジュールは、前述のように、内側ループの繰り返しのたびに、多くの個別の伝達ステップを実行する。図３３は、図３１と３２に示すデータ交換ステップが内部ループの反復回数だけ続けて繰り返される様子を示している。図３４は、図２２のステップ２２０９で、ｐモジュール２８１０で実行する共役ベクトルの計算に関連するデータ交換を示している。ｐモジュールは、スーパー ３４０２から状態ベクトル、制御ベクトルおよび他の計算値を受け取り、これらの値またはこれらの値の一部をｐエージェントに配信し（３４０４〜３４０７）、中間結果をｐエージェントから受け取り（３４０８〜３４１１）、τの値０からＴ_ｋまでの範囲の共役ベクトルを計算してスーパーに返す（値３４１２）。前述のように、現在の好ましい実施形態では、ｐモジュールは共役ベクトルを積分せずに補間するので、図３４に示すような方法で問題を分散する必要はない。

以下の説明では、表記上の都合により、ベクトルｕとｒを併せてベクトルｕとする。つまり、後続の説明では、束縛変数の係数とバリア関数の係数が統合され、唯一の制御ベクトルｕを構成する。本発明の実施形態で採用する反復最適化の戦略には、

の積分が含まれることを想起されたい。ここで、
ｗ(ｙ(τ)，ｕ(τ))＝(Ｈ(ｙ(τ)，ｕ(τ)))^−１Ｇ(ｙ(τ)，ｕ(τ))
ｐ(ｙ(τ)，ｕ(τ))＝(Ｇ(ｙ(τ)，ｕ(τ)))^Ｔｗ(ｙ(τ)，ｕ(τ))
である。一実施形態において、ｄｙ(τ)／ｄτの積分は、順方向オイラー積分法（forward Euler integration method）によって可変ステップサイズΔτ_ｋで、以下のように実行される。

したがって、上の式で明らかなように、ｙモジュールは順方向積分のステップごとにヘッシアン行列の逆行列を求め、ヘッシアンの逆行列とベクトルＧの積としてｗベクトルを計算する必要がある。図３５は、この計算を図式的に示している。ヘッシアン行列の逆行列は、最適化に必要なタスクの中で計算量が群を抜いて最大である。したがって、ｙモジュールはヘッシアン行列の逆行列の計算をｙモジュール・エージェントに分散する。ただし、この分散処理はいくつかの仮定に基づいている。

以下の式を想起されたい。

Ｈ(ｙ(τ_ｋ)，ｕ(τ_ｋ))ｗ(ｙ(τ_ｋ)，ｕ(τ_ｋ))＝Ｇ(ｙ(τ_ｋ)，ｕ(τ_ｋ))

図３６は、ヘッシアン行列とｗベクトルの乗算が分散方式でｇベクトルを生成する様子を示している。図３６に示すように、Ｇベクトルとｗベクトル、およびヘッシアン行列は区分ｉ＝０から５まで縦方向に分割できる。Ｇベクトル３６０２の最初の区分をＧ^０とし、ヘッシアン行列３６０４の最初の縦区分をＨ^０とし、さらにｗベクトル３６０６の最初の区分をｗ^０とする。ヘッシアン行列は、さらに横方向に横区分ｊ＝０から５に分割できる。図３６において、ヘッシアン行列３６０８の最初の横区分はヘッシアン行列内の最初の縦列として表示されている。

ヘッシアン行列を縦と横の両方に区分化すると、ヘッシアン行列がブロック状の部分行列に分割される。たとえば、図３６で、ヘッシアン行列が左上の部分行列Ｈ^０，０（３６１０）から右下の部分行列Ｈ^５，５（３６１２）までの３６個の部分行列に区分化されることを示している。対角線上のブロック部分行列Ｈ^０，０，Ｈ^１，１，Ｈ^２，２，Ｈ^３，３，Ｈ^４，４、およびＨ^５，５が他の対角線上にない部分行列の要素よりもはるかに絶対値の大きな要素を持つ対角優位（diagonally dominant）ヘッシアン行列の場合は、ヘッシアン行列とｗベクトルとの行列乗算をｗベクトルの区分の行列とヘッシアン行列の対応する縦区分内にある対角線ブロック部分行列との乗算で近似してもよい。たとえば、行列全体の乗算を行う場合は、Ｇベクトルの区分Ｇ^０はｗベクトルの区分とヘッシアン行列の最初の縦区分の対応するブロックとの行列乗算による積の和に等しい。これは、次のように表される。

ヘッシアン行列が対角優位である場合は、対角上にない部分行列は要素の絶対値が相対的に小さいため、ｗベクトル区分の計算における重要度は低い。この状況は、ヘッシアン行列が各区分と他すべての区分との間にある依存関係をできるだけ多く削除する方法で区分化された場合に発生する。したがって、特定のｙモジュール・エージェントは、対応するヘッシアン行列の対角ブロック部分行列の逆行列３８０８を計算し、大きなカッコでくくられた項（図３８の３８１２）をブロック対角線ヘッシアン部分行列の逆行列３８０８にかけることによって、ｗベクトル区分３８１０の現在の値を近似して計算できる。ただし、残りのｗベクトル区分３８１４〜３８１７については前に計算された値を使用する。

図３９は第２のｙモジュール・エージェントによるｗベクトル区分ｗ^１の計算を説明する。第２のｙモジュール・エージェント３９０２は、現在のヘッシアン行列の対角ブロック部分行列Ｈ^１，１と現在のＧベクトル区分Ｇ^１および最新に計算されたｗベクトル区分ｗ^０、ｗ^２、ｗ^３、ｗ^４およびｗ^５をｙモジュール２８０８から取得し、現在のヘッシアンブロック対角部分行列Ｈ^１，１の逆行列を計算してから、計算されたヘッシアン行列の逆行列を使用してｗベクトル区分ｗ^１の現在の近似値を計算する。次の計算ステップで、ｙモジュール・エージェント３９０２は、ｗベクトル区分の現在の近似値を他のｙモジュール・エージェントに配信できるようにｙモジュール２８０８に提供でき、次にτが増分されたときに次のｗベクトル区分ｗ^１の値を計算できるように、他のｙモジュール・エージェントが最新に計算した他のｗベクトルの区分ｗ^０、ｗ^２、ｗ^３、ｗ^４およびｗ^５の値を受け取ることができる。この方法を使用すると、ｙモジュール・エージェントは各ステップで同期を取りながらｗベクトル区分の近似値を計算できる。

次に、スーパーモジュールの制御流れ図について説明する。図４０はスーパーモジュールの最も高い概念レベルを示す制御流れ図である。ステップ４００２で、スーパーモジュールは、自らが親軌道から生成された複製または親軌道から分岐した子軌道であるかどうかを判断する。そうでない場合、したがって現在の軌道がルート軌道である場合は、スーパーモジュールはステップ４００４で最適化の初期点を決定する。ステップ４００６で、状態ベクトル、共役ベクトル、および他の値のセットを含むさまざまな計算エンティティを初期化し、さらにブール変数「critical」および「precritical」のブール値ＦＡＬＳＥへの初期化を含む多くの変数の初期化も行う。ルート軌道による初期化に続いて、あるいは複製された軌道の場合にステップ４００２からの制御の流れに従って進んだ結果として、スーパーモジュールは次にステップ４００８で限界Ｔを選択する。さらに、スーパーモジュールは前述の最適化法における連続する内側ループの繰返しを実行する関数「ｕおよびｙをＴまで積分（integrate_u_and_y_to_T）」を呼び出す（ステップ４０１０）。次にステップ４０１２に進み、スーパーモジュールはステップ４０１０のルーチン「ｕおよびｙをＴまで積分」で実行した内側ループの繰り返しに基づいて決定を行うために、ルーチン「評価（evaluate）」を呼び出す。ルーチン「評価」は、前述の最適化法の外側ループの繰り返しをさらにもう１度実行するかどうかを決定するために、ステップ４０１４でスーパーモジュールが評価するブール値「継続（continue）」を返す。実行する場合は、スーパーモジュールは、計算した共役ベクトルの値を返すためにルーチン「ｐを積分（integrate_p）」４０１６を呼び出してから、ステップ４００８に戻って新しい外側ループの繰り返しを開始する。スーパーモジュールが外側ループの繰り返しを継続しないと決定された場合は、スーパーモジュールはステップ４０１８でルーチン「報告（report）」を呼び出してから復帰する。

次に、ステップ４００４、４００８、４０１０、４０１２および４０１８を示す制御流れ図について説明する。ステップ４０１６で呼び出されるルーチン「ｐを積分」と、ルーチン「ｕおよびｙをＴまで積分」から呼び出されるルーチン「ｕを積分（integrate_u）」および「ｙを積分（integrate_y）」は、明確かつ簡潔にするために、制御流れ図を使用せずに数理的に記述する。

図４１は、図４０のステップ４００４で呼び出されるルーチン「初期点（initial_point）」の制御流れ図である。初期点ルーチンは、状態ベクトルｙ_０の初期値を決定するために、ステップ４１０２でユーザーから、あるいはオプティマイザの数理モデル作成からまたは前処理の部分から提供されるすべてのヒントまたは経験則を使用して値を推測する。次に、初期点ルーチンは、ステップ４１０４〜４１０７で構成されるｆｏｒループの中で、各制約ｇ_ｉを評価し、最初に推測した状態ベクトルの点ｙ_０で評価された制約が０より大きいかどうかを判定する。大きい場合は、ステップ４１０７で区分Ｉに制約が指定される。大きくない場合は、区分Ｊにステップ４１０６で制約が指定される。さらにステップ４１１０で、初期点ルーチンは目的関数を次のように構成する。

これは、ステップ４１１２でオプティマイザを再帰的に呼び出し（図３０）、修正された初期状態ベクトルの値を生成することによって最小化される。目的関数Ｆ（上部にチルダ）は、初めは状態ベクトルの初期値を包み込むために目的関数Ｆの問題領域境界を拡張すると考えることができ、次にＦ（チルダ）を最適化することによって問題領域境界を縮小し、初期点を強制的に目的関数Ｆの問題領域の内部に含める。ステップ４１１４で、初期点ルーチンは変更した状態ベクトル値が実現可能かどうかを判断する。実現可能な状態ベクトル値は、局所的最小値を検出する最適化法で局所的最小値に十分近づく可能性の高い問題領域内の点である。変更した状態ベクトル値が実現可能な場合は、最適化の初期点として返される。それ以外の場合は、制御フローはステップ４１０２に戻り、状態ベクトルの初期値を決定する他のパス（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）のステップに関してΔτを決定する。Δτの計算が成功すると（ステップ４６２６で決定する）、ルーチン「ｙを積分」は前進積分ステップを実行し（ステップ４６２８）、必要に応じて束縛変数の係数を更新して（ステップ４６３０）復帰する。Δτの計算が失敗すると（ステップ４６２６で決定する）、ルーチン「ｙを積分」はエラーで復帰する。

オイラー順方向積分の技術について前述したが、形式的には次のように表される。ρの広域的な値は、現在の内部ループ・ステップでｙモジュール・エージェントによって計算された区分化されたρの値から、次のように組み立てられる。

次に、ｙモジュール・エージェントによってｗベクトルの区分の値が収集され、組み立てられてｗベクトルの値が生成される。さらに、後述する方法で連続型反復の増分Δτが計算される。最後に、状態ベクトルは現在の連続型反復の値τ_ｋから次の連続型反復の値τ_ｋ＋Δτに次のようにして伝搬される。

最終的に、連続型反復の値τ_ｋは次のように更新される。
τ_ｋ＝τ_ｋ＋Δτ_ｋ
ただし、ｙモジュール・エージェントの局所的なρの絶対値であるρ^ｉが特定の定数パラメータを下回る場合は、ｙモジュール・エージェントは、ｙモジュールのベクトルＧ区分に関して臨界的とみなされることに留意されたい。ｗベクトルの区分ｗ^ｉとｙモジュール・エージェントの局所的なρであるρ^ｉの計算については以下で説明する。

図４７は、図４６のステップ４６２５で呼び出すルーチン「Δτを計算（compute Δτ）」を示す制御流れ図である。ステップ４７０２で、Δτの見積もりまたは推定値が決定する。この見積もりまたは推定値は、定数パラメータから求めることができる。または、より複雑な実施形態では、以前に計算したΔτの値と状態ベクトルの値、制御ベクトルの値、および他の同様の情報の最新の履歴から、経験的な方法（empirical methods）または経験則（heuristics）によって導くこともできる。説明されている実施形態において、オイラー順方向積分のステップで制御ベクトルがほぼ一定の値をとる場合は、次のようにＦがステップごとに１ずつ減少する。

ｙモジュール・エージェントのいくつかがモースモードの場合は、ｄＦ／ｄτが値「−１」を外れてもよい。しかし、ｄＦ／ｄτの値が０を上回る場合は、ルーチン「Δτを計算」のステップ４７０４で検出するように、エラー状態が発生する。これは、ｄＦ(τ_ｋ)／ｄτが降下の方向ではないので、オイラー順方向積分では実現可能なｙ^＊（τ）の軌跡を生成できないためである。それ以外の場合は、Δτの推定値Δτ'がオイラー順方向積分のステップ式、

に挿入され、状態ベクトルがｙ（τ_ｋ＋Δτ'）に伝搬される。この計算された次の状態ベクトルの値についてその実現可能性がチェックされる（ステップ４７０６）。状態ベクトルの値が実現可能でない場合は、Δτの推定値Δτ'が減少し（ステップ４７０８）、新しい前進状態ベクトルの値が再計算され、その実現可能性がチェックされる（ステップ４７０６）。実現可能な前進状態ベクトルと対応するΔτの推定値Δτ'が検出された場合は、ステップ４７１０でArmigoの条件、

がΔτの推定値Δτ'に関連して評価される。Ａｒｍｉｇｏの条件が満足された場合は、Δτの推定値Δτ'が新しく計算されたΔτとして返される。それ以外の場合は、Δτが近似される。複数の異なる近似方法を利用できる。Δτの以前の実現可能な推定値Δτ^Ａがある場合は、次のようにτｋを中心とする目的関数Ｆの２次近似を作成できる。

次に、Δτの推定値が次のようにして得られる。

Ｆを最小化する降下の方向に沿ったΔτは、Ｆ（チルダ）を最小化するΔτを検出することによって近似できる。これは次の場合に発生する。

したがって、次の式が得られる。

Δτ’＝―ｂ／２ａ

現在のτ_ｋで、以前に計算されたΔτの２つの実現可能な推定値が既知の場合は、３次近似法を利用してもよい。ただし、ステップ４７１２でΔτが近似され、新しい近似値Δτを使用して、状態ベクトルの前進した値が再計算されると、ステップ４７０６でその実現可能性がチェックされる。

前述のように、ｙモジュール・エージェントが臨界点に近く、したがってモースモードの場合は、最適な最適化の軌跡ｙ^＊を検出するときに、前述の標準的なオイラー順方向積分の手続きは役に立たない。一般に、臨界点の近傍では、ｙモジュール・エージェントｉに対応するヘッシアンのブロック対角（block-diagonal）部分行列Ｈ^ｉ，ｉは固有値が負であり、したがって、計算される方向は降下の方向ではない場合もある。しかし、ヘッシアン部分行列Ｈ^ｉ，ｉは、正定値ヘッシアン部分行列

を生成し、状態ベクトルのモースモードの前進積分を実行するように変更できる。

状態ベクトルの前進積分へのモースモードのアプローチの派生は、３つの仮定、すなわち（１）目的関数Ｆ（ｙ（τ），ｕ（τ））は、ｕ（τ）に関してほぼ一定である、（２）ヘッシアン行列Ｈ（ｙ（τ），ｕ（τ））は対角中心であり、これはｙモジュール・エージェントが区分化され、ある程度独立していることを示している、（３）状態ベクトルｙ^ｉ（τ）は、ｙモジュール・エージェントｉがモースモードの間はほとんど変化しない、に基づいている。標準的な方法で計算されたρ^ｉ（ｙ（τ），ｕ（τ））の値で決定されるように、

で、ｙモジュール・エージェントｉがモースモードすなわち局所的な臨界状態に入ると仮定する。前述の第１と第３の仮定に従って、領域内の点

に関して、テーラー近似式（Taylor approximation）により、目的関数Ｆは以下に示すように近似できる。

３次の項を切り捨てると、目的関数は次のように近似できる。

これで、領域内の点

を中心とするブロック対角のヘッシアン部分行列Ｈ^ｉ，ｉは、弱いモース変換を使用して次のように正定値行列に変換できる。

まず、図４８に示すように、ヘッシアン行列は次のようにスペクトル分解によって分解できる。

ＮｘＮの行列Ρ ４８０２は、ヘッシアン行列の固有値を列として備えており、ＮｘＮの対角行列Λ ４８０４には対応するヘッシアン行列の固有値が含まれる。

ヘッシアンを変換する１つのアプローチは、行列Ρの代わりに新しい行列

を作成することである。ただし、この行列の列ベクトルは次のように定義される。

さらに、対角エントリ

を用いて、新たな行列

が定義される。

以前の結果を用いると、正の定符号（positive definite）である変更されたヘッシアン行列は、変更したＰ行列および変更したΛ行列から次のようにして得られる。

さらに、変更されたヘッシアン行列の逆行列は、次のようにして計算できる。

逆行列を必要とするのは、変換されたヘッシアンが特異の可能性があるためである。

ヘッシアンを変換する第２のアプローチは、強モース変換（Strong-Morse Transform）として知られており、モースモードの順方向積分ステップで使用される。この方法では、Λを変更した新しい行列が次のようにして生成される。

そして、変更されたヘッシアンは以下のように算出される。

この変換されたヘッシアンは正則（non-singular）であり、したがって、変換されたヘッシアンの逆行列は次のように明確に定義される。

代替の実施形態では、弱モース変換（weak Morse transform）をモースモードでない点、すなわち臨界点から遠い点に使用できる。

変換されたヘッシアン、

を使用すると、目的関数の近似は次のようになる。

状態ベクトルｙ^ｉ（τ）の第ｉ区分に関してこの近似目的関数の導関数をとると、次のようにベクトルＧの第ｉ区分が生成される。

さらに、τ_ｋで評価されたベクトルＧの第ｉ区分は次のようになる。

したがって、ヘッシアン部分行列は次のように指定される。

繰り返すが、非対角ヘッシアン部分行列は強モース変換では変更されない。ヘッシアンは対角中心であり、非対角部分行列は目的関数の近似にはほとんど影響しないためである。前述の仮定２を使用すると、τ_ｋにおけるブロック対角ヘッシアン部分行列は、τ_ｉ（上部にハット）におけるブロック対角ヘッシアン部分行列にほぼ等しい。

同様に、第３の仮定は以下を意味している。

したがって、τ_ｋで評価された近似目的関数

のベクトルＧは、次のように適切に近似される。

で評価された近似目的関数のヘッシアン部分行列は、次のように近似される。

このような近似によって、エージェントｉに関するｗベクトル区分と局所的なρの近似計算ｗ^ｉとρ^ｉが次のようにして得られる。

図４６のステップ４６３０で、束縛変数の係数を調整してもよい。図１９に示すように、束縛変数の係数は連続的に調整されるが、本発明の少なくとも１つの実施形態において、調整は２進であり、以下の場合に発生する。

ただし、ｙ_Ｎはエンベロープ変数ｚである。
束縛変数の係数が調整される時点では、当初の小さな定数値から特定の大きな定数に変更される。この点で、状態ベクトルも調整する必要があるので、束縛変数の係数を調整しても、状態ベクトルを実現不可能な点に強制的に配置することはない。状態ベクトルは次のように調整される。

ただし、パラメータ「ipm_c_offset2」の値は１に近い。

ｕ＿ａｎｄ＿ｒモジュールは、内部ループの繰り返しごとに制御ベクトルｕの次の値を計算するのに必要である。この計算は、ｙモジュール・エージェントにおけるｙベクトルの計算によく似ている。前述のように、次の制御ベクトルｕを計算するための繰り返しの式は以下のようになる。
ｕ_ｋ＋ｌ＝ｕ_ｋ−φ_ｋ

静的な問題を定式化するために、ここでは繰り返しごとにハミルトニアンの凸２次近似（convex quadratic approximation）を使用する。ここでは最小値を表す明示的な式を検出できるので、必ずしもパラメータσを使用する必要はない。ｕを計算するには、唯一のステップが必要である。

ここで、ｕ_ｋ＝ｕ|_τｋを、関数Ｈ（ｙ，ｕ，ｐ）の点とする。ｙ_ｋとｐ_ｋが固定された場合に、二次テーラー展開を使用したｕ_ｋの周辺の二次近似は、次のようにして計算できる。

ここで、Ψは∂Ｈ／∂ｕとして定義される。これは２次関数であるが、必ずしも凸ではない。∂Ψ／∂ｕが必ずしも正定値（positive definite）ではないためである。モース変換を使用して、∂Ψ／∂ｕを正定値の凸近似

に変換できる。これは、前述のヘッシアン行列の変換と同様である。

これで、ハミルトニアン行列を次のように近似できる。

Δｕ＝ｕ−ｕ_ｋと定義すると、次のようになる。

以下の説明では、ｙ_ｋ、ｕ_ｋ、ｐ_ｋへの依存性については明確に示していない。上の関数の最小値は、目的関数の勾配を０に設定し、Δｕについて解くことによって検出できる。したがって、

となり、次の式が得られる。

ここで、上付き文字Ｒは前述のように反転を表す。ここで、

と定義すると、次の式が得られる。

これで、次の制御演算の反復更新が可能になる。
ｕ_ｋ＋１＝ｕ_ｋ＋Δｕ
しかし、Ψおよび∂Ψ／∂ｕについて解くのは複雑な計算プロセスである。ハミルトニアンは次のように記述される。

チェインルール（chain rule）と、制御演算ではｐが定数として扱われるという事実とを使用すると、α＝０，．．．，Ｍ−ｌにおいて、

となり、α，β＝０，．．．，Ｍ−１において、

となる。

ここで、α＝０，．．．，Ｎ−１において、次の定義、

により、このｕに関する１次導関数は次のように表される。

ただし、α＝０，．．．，Ｎ−１、β＝０，．．．，Ｍ−１である。

ヘッシアン行列の逆行列Ｈ^−１は次のように定義される。
Ｈ^−１Ｈ＝Ｉ
さらに、行列の要素について次のように記述できる。

ただし、β＝０，．．．，Ｎ−１であり、さらに

である。両辺で積の法則を使用してｕに関する導関数をとると、次の式が得られる。

ただし、α，η＝０，．．．，Ｎ−１、γ＝０，．．．，Ｍ−１である。移項すると、次のようになる。

両辺にＨ^−１をかけて次の式を得る。

式を簡素化すると、次のようになる。

Ｈ^−１の要素を表す明示的な式を次に示す。

Ｈ^−１の式を

に代入すると次の式が得られる。

ただし、α＝０，．．．，Ｎ−１、β＝０，．．．，Ｍ−１である。

２次導関数は次のようになる。

Ｆの勾配とヘッシアン行列はｕに関して線形なので、これらの項のｕに関する２次導関数は０であるという事実によって式が簡素化されることに留意されたい。Ｈ^−１の式を上の式に代入すると次の式が得られる。

これは次のように簡素化できる。

次に、分母（denominator）ρが次のように定義される。

ｕに関する１次導関数と２次導関数は、積の法則を使用して次のように表される。

ただし、α＝０，．．．，Ｍ−１、

ただし、α，β＝０，．．．，Ｍ−１。

次に、形式的な分散アルゴリズムについて説明する。

の両辺に

をかけて次の式を得る。

ｒ番目のエージェントの項における関係は次のように表される。

ただし、ｉ＝０，．．．，Ｍ−１。移項して次の式を得る。

以下の表記

を使用して、両辺に

の逆行列をかけると次の式が得られる。

したがって、次の式が得られる。

したがって、制御を次のように更新できる。
ｕ^ｒ _ｋ＋１＝ｕ^ｒ _ｋ＋Δｕ^ｒ 。

ここで、エージェントの項に以下に示すハミルトニアンの定義を記述する。

これで次の式が得られる。

ここで、

（ただし、α＝０，．．．，ｎ^ｉ）、および

を考慮すると、エージェントＩのｗベクトルの各要素は次のように表すことができる。

前の式にＨ^−１の定義を代入すると、次のようになる。

移項すると次の式が得られる。

ＧとＨのｕに関する２次導関数は０なので、次のように表される。

ただし、

ここで、項を分割し、終わりの２項の添字を変更すると、次のようになる。

したがって、次の式が得られる。

さらに、ＧとＨのｕに関する２次導関数は０であるという事実を使用して、次の式を得る。

最終的に次の式が得られる。

さらに、この導関数は次のようになる。

共役変数の伝搬は、状態変数と制御変数の伝搬とは異なっている。状態変数と制御変数は、密結合し、共に前方に伝搬する。対照的に、共役変数は状態変数と制御変数の結合とは関係なく、適切な時期に後方に伝搬する。このためには、ｙモジュールとｕ＿ａｎｄ＿ｒモジュールで計算した値のデータベースが必要である。ｐの微分方程式はｐに関して線形なので、解は積分でなく補間法によって近似できる。多くのアプリケーションで、ベクトルｐの値は収束限界周辺の過渡現象（transient around the convergence horizon）を除いてほとんど一定である。ｐの式における行列Ａもｐに依存しない。

ｐを計算するために、ｙモジュールとｕ＿ａｎｄ＿ｒモジュールで計算されたｙとｕの軌跡を表す一連のノードの値が作成され、これをｐモジュールで使用する。まず終点の条件ｐ（Ｔ）が計算され、次にｙの初期の計算まで後方に補間するループの中で、ｐの初期の値が計算される。

実用的にするために、ここでは行列Ａがほとんど一定であると仮定される場合は、τ_ｋ−１−τ_ｋ（ただし、τ_ｋ−１＜τ_ｋ）のサイズの増分に問題を分割する。次に、後退積分を実行すると、τ_ｋが指定された場合にτ_ｋ−１における解は次のように表される。

ｐ(τ_ｋ−１)＝exp(Ａ(τ_ｋ−１−τ_ｋ))ｐ(τ_ｋ)

この式の指数は、以下の近似によって解決される。

ただし、Ｎは近似の次数である。

修正（Repair）
ここまで、一般的な最適化の特定の計算方法について説明してきた。しかし、現実世界の最適化における付加的な側面は、実際の時間区間にわたって繰り返される最適化に関連する。たとえば、メーカーなどの企業では、最適化の技術を使用して、従業員のスケジューリング、パーツや原料の注文、パ―ツ、原料、中間製品の製造工場内での配布、そして特に顧客への製品の配送を含む製造工場の運転に関するさまざまな側面をすべて最適化しようとする。一般に、こうした企業の最適化を行うための非常に高次元の問題領域では、製造工場運転中に発生する管理されない多くのイベントや変化によって、時間が経つと共に変化する条件に合わせて製造工場の制御を再調整するために、最適化を何度も実行しなければならない場合がある。最適化以外に、同様の問題は他の多くの計算タスクにも発生する。たとえば、複雑な系のモデル化では、系に関してランダムに発生する様々なイベントや混乱に対応するために、系の状態に関する中間的な計算を導入し、変化する環境への系の実際の時間応答をそのモデルが反映するようにしなければならない場合がある。

図４９は、こうした反復最適化の戦略を示している。図４９と後続の図５０および５２〜５６において、横軸４９０２は時間に対応する。図４９では、状態変数と制御変数を含む系の準最適な状態関数が、図４９に縦の破線で示す各時点（たとえば縦の破線４９０４）で、最適化および／またはモデル化の技術を使用して計算される。計算された状態関数は、図４９には曲線の一部分として表されている。たとえば、曲線の一部分４９０６は、時刻４９０８の点で開始されており、この時点で状態関数の計算が開始され、別の時刻４９１０の点まで展開されている。もちろん、計算された状態関数は、現実の世界に適用する上で、超次元空間における複雑な制御の軌跡であるが、説明を明確にするために、曲線の一部分４９０６のような２次元の制御の軌跡を使用している。斜めの破線（たとえば斜めの破線４９１２）は、前述の連続型反復変数τに関する積分を使用した状態関数の計算を表しており、状態関数を計算する最適化技術の計算限界Ｔに関連する特定の長さの時間４９１４を必要とする。新しい状態関数、たとえば、時刻４９１８で開始され、時刻４９２０で終了する時間の範囲に適用できる新しい状態関数４９１６が計算されると、この新しい状態関数は直前の状態関数（この例では状態関数４９０６）と一部の重複する部分で比較的接近してはいるが、一致してはいない。このように、連続して計算された２つの状態関数の間には、ずれ４９２２または不連続が存在する。この不連続は、前の状態変数が計算された時刻と現在の状態変数が計算された時刻との間に発生した系の条件や環境の変化によって発生する。この場合は、時刻４９０８で開始され、時刻４９１８で終了する時間区間にイベントや変化が発生している。

残念ながらこうした不連続は系の混乱を表している。たとえば、製造工場において、不連続は従業員のスケジュール変更、配送用トラックの出発時刻と行き先の変更、およびその他の同様の変更を表していてもよい。しかし、複雑な系は慣性（inertia）をもつ傾向があるので、大きな混乱を克服するのは困難であり、前回計算した最適な状態関数が最適でなくなる段階まで、自らが系の環境をある程度変更し、複雑なシステムをさらに混乱させるリプルやカスケードの効果を発生するという問題を克服するには、長い時間と多くの資源が必要になる場合がある。

最適化の繰り返しによって発生する不連続を改善する１つの可能なアプローチは、最適化を実行する間隔を短くすることである。図５０は、図４９に示す状態関数計算による最適化の間隔を短縮したものを示す図である。状態変数をより効率的に計算することによって、図５０では状態関数の計算を表す破線の傾きが大きくなっており（たとえば破線５００２）、最適化を短い間隔ΔＴ ５００４で実行してもよい。しかし、その結果として混乱の頻度が増し、混乱の頻度が増したこと自体によってリプルやカスケードの効果が発生し、場合によっては不安定性がますます増大することにもなり得る。一部の事例では、非常に短い時間区間に特定の種類のイベントが同時に発生したことによる不足減衰（under damping）の結果として、より大きなずれが発生することさえある（たとえば、ずれ５００６）。

図５１は、最適化の繰り返しによって発生するもう１つの結果の類推を示している。図４９〜５０では、新しい状態関数は、前に計算された状態関数による制御を終了し、新しく計算する状態関数による系の制御を開始する点を除いて、前に計算された状態関数に関係なく新たに計算されることを想定していることに留意されたい。しかし、前の状態関数計算の履歴を考慮しないと、長い時間が経過した場合に、極めて非最適な状態関数の解が得られる可能性がある。同様に、図５１に示す単純な地形図で表される地域を点Ａ ５１０２から点Ｂ ５１０４まで移動する問題について考察する。点Ａから点Ｂに至る第１のパス５１０６は、ある人がコンパスだけを頼りにたどるであろうパスを表しており、コンパスの針の動きが止まるのを待たずに、瞬間的なコンパスの読みだけに基づいて頻繁に軌道修正を繰り返している。こうしたパスは、点ＡとＢの間の移動に関して、時間的にも距離的にも全くの非最適解を表している。

第２のパス５１０８は、コンパスの動きが落ち着いてから読む必要があることを認識している人、がけに磁鉄鉱が露出している場合などの局所的な条件によって混乱が発生する可能性があることを認識している人、および障害物が現れた場合は、コンパスの示す方向でなくその障害物に慎重に対処することを認識している人がたどるであろうパスを表している。第２のパス５１０８は、点ＡとＢの間の移動に関して、はるかに効率的かつ効果的なパスを表しているのは明らかである。

第３のパス５１１０は、地形数学（topographic math）を理解でき、コンパスを正しく使用でき、川にかかった橋は目的地への重要な中間点であることを認識している人、森の道に沿って歩く方が道を開拓するよりはるかに容易なことを認識している人がたどるパスである。図５１から類推すると、過去の情報に基づくことができる人、単純で限られた現在の状態以外の情報に基づくことができる人、追加の情報に基づいて将来的な解を予測できる人は、点ＡとＢの間の最適なパスを検出する可能性の高い人である。図４９〜５０に示す反復最適化問題に戻り、局所的、瞬間的な状態関数を考慮するのみでなく、前に計算された状態関数を使用してより長い期間にわたってより多くの情報を取得する方法を検出することが望ましい。

図４９〜５０に関連する前述の考慮事項の結果として、長期にわたる履歴情報と最新の予測を使用してある程度連続型最適化のアプローチは、最適化を目的とする本発明の実施形態に追加される機能、「修正に基づく連続型最適化のアプローチ」（ＲＣＯ：Repair-based Continuous Optimization）を表している。図５２は、変動時間枠（sliding window）の概念を示している。図５２に示すように、時刻ｔ_１−ΔＴ ５２０４で始まり、時刻ｔ_１＋Ｔ−ΔＴ ５２０６で終了する時間長Ｔの前の時間枠５２０２は、状態関数を計算するための参照フレームとして使用される。この時点から、時間長ΔＴ ５２０８が経過している。したがって、時間長Ｔの新しい現在の時間枠５２１０は、現在の時刻ｔ_１ ５２１２における参照フレームを提供する。

図５３に示すように、前の時間枠５３０２で計算された状態関数は時刻ｔ_１−ΔＴ ５２０４からＴ_１−ΔＴ＋Ｔ ５２０６に展開される。図５４は、前に計算された状態関数の現在の時間枠への展開を示している。図５４に示すように、前に計算された状態関数５３０３は、５４０２から現在の時間枠の終わり（時刻ｔ_１＋ΔＴ ５４０４）まで展開される。図５５は、現在時点における考慮事項、過去のイベントに関する考慮事項、および現時点で前に計算された状態関数が陳腐化し、最適でなくなる可能性のある環境の変化、および現在のウィンドウで表される現在の参照フレームで前の状態関数の情報を新しい最適な状態関数の検出に組み込む方法を示している。

図５５において、縦の線５５０２〜５５０５は、状態関数の計算時間中、すなわちｔ_１−ΔＴおよび現在時刻ｔ_１に発生したイベントを表している。こうしたイベントによって、系の状態や環境が変化し、前に計算された状態関数が最適でなくなる可能性がある。たとえば、製造工場で、ストライキやインフルエンザ・ウィルスによって多くの従業員が製造工場内の作業を実行するようにスケジューリングできなくなる可能性がある。したがって、前に計算された最適なスケジューリングは、従業員不足に直面した場合は完全な次善最適解であろう。図５５に示すように、状態関数の計算の基準となる時間依存パラメータに関するさまざまな予測は、時刻ｔ_１では変化している可能性がある。こうした予測の変化は、破線５５０６で表されている。たとえば、製造工場のコンテクストにおいて、主要な顧客がこの顧客が通常購入する数の２倍の製品を緊急に要求してきたり、競合他社で製造された新しい製品の出現によってこの製造工場で製造する製品の需要が突然減少したりする可能性がある。このように、前に予測した製品需要が現在および将来の需要を反映しない場合がある。

さらに、状態関数の基準となるさまざまな考慮事項やパラメータに対して当初指定された優先順位が変化する可能性がある。図５５の点線５５０８は、現在の参照フレームで発生した優先順位の変化を前の参照フレームに関連して示している。たとえば、木曜日に状態関数が計算された時点で、時刻ｔ_１−ΔＴからｔ_１までの間に、差し迫った道路の閉鎖に関する情報が発表されていなかったとする。木曜日の道路の閉鎖に関する情報によって、南側に住む従業員が通勤時に橋の閉鎖の影響をほとんど受けない時間帯で作業できるスケジューリングの優先順位が大幅に向上する。

さらに、系の動作における重大な混乱を回避するために、系は新しく計算された状態関数が特定の区間にわたって同じ状態を示しているかどうかにかかわらず、直前に計算された系の状態に従う必要がある。たとえば図５５で、それぞれＴ_２ ５５１４とＴ_４ ５５１６で始まり、それぞれＴ_３ ５５１８とＴ_５ ５５２０で終了する２つの時間枠５５１０と５５１２は、新しく計算された状態関数内に状態関数の古い値５５２２と５５２４を保持する必要がある。このように、前の状態関数の計算で得られる履歴の情報は、現在の参照フレームに保持される。このようなフェンス区間は、前に計算された状態関数に基づいて確立するシステム制御における大きな障害や混乱を避けるために必要であろう。

イベントに関する情報、環境や優先順位の変化に関する情報、および前に計算された特定の状態を保持する必要があるかどうかを使用して、新しい状態関数を計算せずに、前に計算され、展開された状態関数を変更し、現在の時間枠に関する準最適解を生成することができる。図５６は、現在の時間枠に関する準最適な状態関数を得るための、前に計算された状態関数の修正を示す図である。図５６に示すように、前に計算された状態関数の修正には、現在の時間枠内の時刻における、前に計算された状態関数と現在の時間枠に関する新しい準最適な修正された状態関数との差分の計算が含まれる（たとえば差分５６０２）。差分５６０２のような差分の計算では、前の時間区間ｔ_１−ΔＴ ５２０８に発生したイベント、環境や優先順位の変化（図５５の５５０８と５５０６）、特定の時間区間５５１０と５５１２に関して前に計算された状態５５２２と５５２４を保持する必要があるかどうかが考慮される。したがって、修正に基づく連続型最適化技術によって、系の制御における混乱や障害を最小化し、不連続な制御を回避するように最適化全体を繰り返し実行することなく、ある程度連続型状態関数の再計算が可能になる。

以下は、本発明の１つの実施形態を表す修正方法の、より形式的な記述である。制御最適化問題は、その解が前述の修正方法で修正される制御状態関数を提供し、以下のように指定される。

ただし、制約

を受け、ここで、ＶはＲ^ｒのコンパクト・サブセット（compact subset）であり、ｄ（ｔ）は駆動（driving）関数（たとえば、決定論（deterministic）または予測需要（forecasted demand））である。

および｛ｇ_ｊ，ｊ＝１，．．．，ｍ｝は、十分になだらかである。以下のようなスカラー変数ｘ_ｎ（ｔ）≧０が導入される。

これで、最適制御問題は次のようになる。

ただし、ｔ∈［ｔ_１，ｔ_１＋Ｔ］の場合

を受け、ここで、

である。こうした問題の解ｖ^＊（ｔ）は、「最適制御（optimal control）」と呼ばれており、対応する軌跡ｘ^＊（ｔ）は最適な軌跡（optimal trajectory）と呼ばれている。前述の制御状態関数（たとえば、図４９の４９０６や４９１６）は、ｘ（ｔ）とｖ（ｔ）（ｔ∈［ｔ_１，Ｔ］）との組合せである。

あらかじめ計算された制御状態関数は、図５２〜５６で使用して表すと、次のｖ^{（ｔ１−ΔＴ）＊}（ｔ）およびｘ^{（ｔ１−ΔＴ）＊}（ｔ）の組合せとなる。計算された差分（たとえば、図５６の５６０２）は、「修正（the repair）」と呼ばれ、δｘ^ｔ１ _ｌ（ｔ）およびδｖ^ｔ１ _ｋ（ｔ）の組合せとして定義される。フェンス・インディケータ（Fencing indicator）関数は、制御状態変数のそれぞれについて、ｔが修正された制御状態関数に偏差が許容されない区間内にあるかどうかを指定する。

時間区間［ｔ_１，ｔ_１＋Ｔ］内の有限個のフェンス区間（fenced intervals）を想定すると、（ｉ＝１，．．．，ｎ）（ただし、ｎはモデルの次元）はモデルの次元であり、増大する状態の軌跡δｘ^ｔ１ _ｉ（ｔ）は次のように指定される。

ここで、（ｔ_１−Δｔ）^＊の表記は、偏導関数が参照軌跡に沿って評価されることを示している。

集合｛[ｔ_ｓ ^ｉ，ａ、ｔ_ｓ ^ｉ，ｂ]、ｓ＝１,...,Ｍ_ｉ｝は、状態ｉ＝１,...,ｎのフェンスのある時間の部分区間である。時間区間［ｔ_１−ΔＴ，ｔ_１］内に発生するすべてのイベントＡ^ｅ _ｉδ(ｔ−ｔ_ｅ)は、次のようにｔ_１に発生する唯一のイベントＡ^ｔ１ _ｉδ(ｔ−ｔ_１)に近似される。

δｘ^ｔ１ _ｉ（ｔ_１）＝ｘ^ｔ１ _ｉ（ｔ_１）−ｘ^{（ｔ１−ΔＴ）＊}(ｔ)−Ａ^ｔ１ _ｉは、δｘ^ｔ１ _ｉ（ｔ_１）の初期条件を指定する。ただし、ｖ^{（ｔ１−ΔＴ）＊}（ｔ）およびｘ^{（ｔ１−ΔＴ）＊}（ｔ）は、［ｔ_１−ΔＴ，ｔ_１＋Ｔ−ΔＴ］上で定義される。前述のように、制御状態関数は区間［ｔ_１，ｔ_１＋Ｔ］に展開する必要がある。展開は一定の制御によると仮定し、また展開内のフェンスは許可されないと仮定すると、次式のようになる。

ここで、ｔ∈［ｔ_１＋Ｔ−ΔＴ，ｔ_１＋Ｔ］（ｉ＝１,…,ｎ）であり、

である。

区間［ｔ_１，ｔ_１＋Ｔ］に関する動力学（dynamics）は、次のように与えられる。

Φ（ｘ（ｔ_１＋Ｔ））で

を表すと、２次までの展開によって次の式が得られる。

修正の目的は、上の基準を最小化しながら、同時にノミナル（nominal、名目上の、）軌跡に関する変化を最小化することである。この基準の組合せは次のとおりである。

ここで、ＱとＲは企業の基準（enterprise criterion）を満たすことと変化を最小化することとの兼ね合いを定義するためにユーザーが指定する定数の正定値行列である。

修正問題の制約は次のとおりである。

この式を１次まで展開すると次のようになる。

制約の不等式を満たすために、次の式が与えられる。

展開された制約の第２項と第３項は次を満たす必要がある。

メタ制御と階層的プロセスの構築
これまでの記述は、最適化と後続の点での最適化法によって生成された制御状態関数の調整に焦点を合わせている。しかし、本発明には、複雑で階層的な計算システムとモデルを構築するため、さらに簡単か複雑かにかかわらず一般的な計算を制御するための、より一般的で強力な方法が含まれる。前述の計算制御によるこの一般的な方法の例を示し、当初の最適化問題に追加の制約を加えてメタレベルの最短時間制御問題の目的関数Ｊを生成し、さらにメタレベルの問題に関するハミルトニアンを最小化することによって準最適制御ポリシーを生成する。

図５７は、計算プロセスのメタ制御の１つのレベル示している。図５７では、最初の計算問題または計算モデルである５７０２は、メタレベル軸５７０４の最も低いレベルに示されている。計算問題、またはモデルは離散的であり、反復

によって形式的に説明できる。前述のように、離散的な計算問題またはモデルに連続型反復変数のパラメータ化を導入することによって、当初の計算問題を連続的とみなすことができ、問題を解決するために強力で連続型計算方法を適用できる。したがって、図５７に示すように、この問題は面５７０６では離散的に思われるが、別の面５７０８では連続型（continuous）ｄｘ（τ）／ｄ（τ）＝ｆ（ｘ（τ）、ｕ（τ））とみなされる。離散から連続への変換は、メタ制御の方法を有効にする重要な機能である。問題の連続型ビューを使用することにより、追加の制約を導入でき、前述の最適化問題で構築された最短時間制御問題などの新しい問題を作成できる。この制約を追加する機能の重要な側面は、新しく追加された制約が限定的ではないこと、または当初の問題への関連付けが必要ないことである。また、前述の実施形態において、初期の問題には最適化問題に対する最適制御解の検出が含まれるが、追加の制約は当初の最適化問題において時間的に効率のよい計算に向けて計算を操作することであり、準最適制御ポリシーの検出に関連する。その他の種類の制約を追加してもよい。たとえば、最適化を実行する曜日を操作するメタレベルの制約を追加してもよい、時間とメモリの両方を効率的に使用して計算するための制約を追加してもよい。このように、追加の制約は当初の問題に重ね合わせるメタレベルのビューを表している。新しい制約を追加できることによって、たとえば、周知のトップダウン型アプローチで問題を解決するといった複雑なソフトウェアの設計で使用される概念的な設計のレベルと同様に、階層的な計算による問題解決方法を容易に適用できる。

図５７では、大きな矢印５７１０で表される変換の最初の部分として、制約の追加が示されている。当初の問題全体が、新しい問題における制約の内部と新しい目的関数の項内の両方で、新しい問題に伝えられていることに留意されたい。したがって、新しい問題は当初の問題を含んでおり、当初の問題の解は新しい問題の解の部分集合になっている。したがって、メタレベルではあるが、新しい制約と当初の問題とは全く異なるビューに焦点を合わせることができ、当初の問題は新しい問題として符号化される。変換の第２の要素では、ハミルトニアンが順に使用され、メタ制御ポリシーを計算するための離散型（discrete）の計算方法を決定するための式を提供する。離散型のメタ制御ポリシーの式５７１２は、新しい連続する繰り返し変数σを導入することで、それ自体を連続型の問題５７１４とみなすことができる。

図５８は、複数レベルのメタ制御および／または複雑な計算プロセスの階層構造を示している。図５８に示すように、当初の離散問題５８０２からのメタレベルの変換を繰り返すことで、メタレベルの軸５８１０に沿って任意の数のメタレベル問題が得られる。前述のように、それぞれの新しいメタレベルにより、当初の問題および低レベルのメタレベルを高レベルのメタレベルへと伝搬させながら、新しい制約を追加し、新しい問題を生成できる。例として、現在科学界で大きな関心がよせられている、生きた細胞をモデル化するという計算によるモデル化の問題を考える。細胞は数多くの概念的なレベルとして階層的に見ることができる。すなわち、量子力学および統計力学によって支配される分子レベルから始めて、超分子複合体（super molecular complexes）、巨大分子構造（macromolecular structures）、およびコンパートメント（compartments）、細胞小器官（organelles）、遺伝子制御（genetic control）、環境の影響、および他の無数の概念レベルへと向上し、それぞれは異なる観点でとらえることができ、さまざまな力と原理で操作される。しかし、各レベルは細胞全体に寄与し、分子レベルのイベントの作用や状態は上記の概念レベルに深く影響している。たとえば、特定の染色体内における核酸塩基から別の核酸塩基への自然変換は、遺伝子の変化を引き起こし、次に重大な組織構造または一時的な制御構造の構成要素である遺伝子産物（gene product）の変化を引き起こす。これは、さらに細胞全体の存続に重大な影響を与える可能性がある。さらに、各レベルの問題のさまざまな部分は、本質的に整数、ブール値、または他の型の変数で表現できるが、問題の他の部分は連続変数の項で定義できる。

本発明の方法は、生きた細胞やリアルタイムの分散オペレーティング・システムなどのシステムの複雑な計算モデルを開発するのに適用できる。前述の非実数変数を実数変数に変換し、さらに浮動小数点変数に変換する方法は、さまざまな下位の問題の本質的な表現を、本発明の方法を適用できるすべて浮動小数点からなる問題空間に組み込む方法を提供する。さまざまなメタレベルのビューを定義できることにより、低レベルのビューを伝搬させながら、高レベルのメタレベルで低レベルのビューを明示的に指定する必要なしに、細胞のモデル化の問題から低レベルのビューの複雑性をすべて保持する階層的なモデルの作成が可能になる。

より形式的には、次のように反復アルゴリズムが提供される。

ここで、ｕ_ｋは制御変数、ｘ_ｋは決定変数である。このアルゴリズムは次の式に収束する。

このアルゴリズムの詳細は、関数ｆ(チルダ)とパラメータｕ_ｋでコード化される。

連続化（continualization）の手続きでは、反復はｘ（τ）とｕ（τ）に関する

となるような制御された微分方程式に変換される。変数τは、離散型の反復変数_ｋに対応する連続型反復変数である。この微分方程式は、すべてのｕ（□）∈Ｕ（□）について次の形をとる。
ｄｘ（τ）／ｄτ＝ｆ（ｘ（τ），ｕ（τ））

ここで、Ｕ（□）は実現可能なパラメータ値の集合である。

最適制御問題は次のように表現される。

関数φとΨは、アルゴリズム設計者がメタ目的（meta objectives）を達成するために選択する。

ポントリャーギン（Pontryagin）の最小原理と、

で与えられるハミルトニアン（ここで、ｐ（τ）は共役（costate）である）を使用すると、ｘ^＊（τ），ｕ^＊（τ），∈［０，Ｔ）が最適であるために必要な条件は次のとおりである。

ここで、ハミルトニアンは次のようにｕ^＊（τ）におけるｕの関数としての絶対的な最小値をとる。
Ｈ(ｘ^＊(τ)，ｕ^＊(τ)，ｐ(τ))≦Ｈ(ｘ^＊(τ)，ｕ(τ)，ｐ(τ))，
∀ｕ^＊(τ)∈Ｕ(τ)，τ∈［０，Ｔ)。

収束数列

が、次式により生成される。

ここで、ｘ（τ）＝ｘ^＊（τ）、τ∈［０，Ｔ)のとき、

となる。

前述の収束数列の生成式は、次に示す当初の反復アルゴリズムに類似している。

ここで、

は、

に類似している。

ｓの連続型パラメータ化をσとして前述の連続化手続きを使用すると、次の式が得られる。

ここで、ｖ（τ，σ）は

を連続化した形であり、Ｗ（ｘ（τ），ｖ（τ），ｐ（τ））は、

を連続化した形である。
要約すると次のようになる。

本発明について、特定の実施形態に関連して説明してきたが、本発明をこの実施形態に限定する意図はない。本発明の精神を逸脱しない変更は、当業者には言うまでもない。たとえば、ほとんど制限なく多種多様である。

以下の説明において、本発明の理解を深めるために、説明を目的として特定の専門用語を使用した。ただし、特定の詳細が本発明を実施するために必要でないことは、当業者には言うまでもない。他の例では、基盤となる発明からの不要な逸脱を回避するために、周知の回路とデバイスがブロック図の形で示されている。したがって、本発明の特定の実施形態に関する以上の説明は、例示および説明を目的としており、すべてを記述するものではなく、本発明を開示した形式のみに限定するものでもない。以上の教示に関連して、さまざまな変更や変形が可能なことは明らかである。本実施態様は、本発明の原理の最適な説明を提供するために選択され、説明されており、他の当業者はこれを実際に適用することによって、本発明とさまざまな変形を伴うさまざまな実施形態を、予定される特定の用途に適した形で最適に利用できる。本発明の範囲は、前述の請求項とその均等物によって定義されるものである。

非常に単純な２次元の最適化問題を示す図である。 非常に単純な２次元の最適化問題を示す図である。 本発明の１つの実施形態を表す最適化の方法およびシステムの最高レベルのビューを示す制御流れ図である。 モデル作成と前処理を説明するための仮説的な問題を示す図である。 仮説的な問題のさらなる側面と特性を、図３の描画法を使用して示す図である。 工場から卸売店、および卸売店から顧客への商品の出荷に関連するコストを示す図である。 仮説的な問題の付加的な特性を示す図である。 ２か所の製造工場、３か所の卸売店の候補、５か所の顧客による、コスト、容量、需要の数値に関連する単純な仮説的な問題を示す図である。 図２に示す高レベルの制御流れ図の前処理のステップを示すより詳細な流れ図である。 図７に示す特定の最適化問題を最初に数理的に表現した記号形式のモデルを示す図である。 仮説的な数理モデルから浮動小数点変数のみを使用したモデルへの変換を示す図である。 標準の形式に変換した仮説的な問題の数理モデルを示す図である。 変換された等式制約に対する束縛変数の効果を示す図である。 変換された等式制約に対する束縛変数の効果を示す図である。 変換された等式制約に対する束縛変数の効果を示す図である。 変換された等式制約に対する束縛変数の効果を示す図である。 等式制約を不等式制約に変換した後の仮説的な問題を記号形式で数理的に表現したモデルを示す図である。 束縛変数の効果を図式的に示す図である。 エンベロープ変数ｚを導入した後の仮説的な問題の数理モデルを示す図である。 バリア関数の概念を示す図である。 バリア関数の概念を示す図である。 バリア関数の概念を示す図である。 バリア関数の概念を示す図である。 バリアの制約を追加した後の仮説的な問題の記号形式による数理モデルを示す図である。 初期点（初期点）から局所的最小値に至る最適化の軌跡の方向を示す問題領域内の勾配フィールドの利用を示す図である。 最適化の軌跡を示す図である。 準最適状態ベクトルの軌跡ｙ^＊と準最適制御ポリシーｕ^＊およびｒ^＊に基づく反復最適化法を図式的に示す図である。 反復最適化法全体を図式的に示す付加的な図である。 反復最適化のステップ（図２の２０３）を説明する制御流れ図である。 椀形の問題領域と多くの中央の突起の面境界を示す図である。 図２３と同じ問題領域を示し、臨界点を明示的に示す図である。 臨界点の拡大表示を示す図である。 第１の検索技術を示す図である。 問題領域の第２の検索戦略を示す図である。 最適化スレッドのモジュールとモジュール編成を示す図である。 スーパーモジュールで保守するデータベースの内容を部分的に示す図である。 図２８の描画法を使用して、最適化法の１つの外部ループの繰り返し中に階層構造をなす最適化軌跡のモジュール間で情報が流れる様子を示す図である。 図２８の描画法を使用して、最適化法の１つの外部ループの繰り返し中に階層構造をなす最適化軌跡のモジュール間で情報が流れる様子を示す図である。 内部ループの次のステップ（図２２の２２０５）におけるデータの流れを示す図である。 図３１と３２に示すデータ交換のステップが内部ループの繰り返しによって正常に繰り返される様子を示す図である。 図２２に示すステップ２２０９のｐモジュール２８１０で実行する共役ベクトルの計算に含まれるデータ交換を示す図である。 ｗベクトルの計算を図式的に示す図である。 ヘッシアン行列とｗベクトルの乗算によって得られるＧベクトルを分割された形で示す図である。 ヘッシアン行列とｗベクトルの乗算を分解して示す図である。 ｙモジュール・エージェントによるヘッシアン行列の逆行列の分解を容易にするＧベクトル区分の計算の数理的な操作を示す図である。 第２のｙモジュール・エージェントによるｗベクトル区分ｗ^１の計算を示す図である。 最高の概念的レベルのスーパーモジュールを示す制御流れ図である。 図４０のステップ４００４で呼び出すルーチン「初期点」を示す制御流れ図である。 ルーチン「次のＴを選択」を示す制御流れ図である。 図４０のステップ４０１２で呼び出すルーチン「評価」を示す制御流れ図である。 ルーチン「報告」を示す制御流れ図である。 図４０のステップ４０１０で呼び出すルーチン「ｕおよびｙをＴまで積分」を示す制御流れ図である。 ルーチン「ｙを積分」を示す制御流れ図である。 図４６のステップ４６２５で呼び出すルーチン「ΔＴを計算」を示す制御流れ図である。 ヘッシアン行列がスペクトル分解によって分解できることを示す図である。 制御問題の最適化法による連続型状態関数の計算を示す図である。 図４９に示す状態関数の短縮された区間による計算を示す図である。 繰り返しによる現在の状態ベースの最適化から得られる結果の類推を示す図である。 変動時間枠の概念を示す図である。 前の時間枠で計算されたＴ_１−ΔＴからＴ_Ｔ−ΔＴまでの状態関数を示す図である。 現在の時間枠の時間区間において前に計算された状態関数の展開を示す図である。 現在時点における考慮事項と、過去のイベントに関する考慮事項と、現時点で前に計算された状態関数が陳腐化し、最適でなくなる可能性のある環境の変化、および現在の時間枠で表される現在の参照フレームで前の状態関数の情報を新しい最適な状態関数の検出に組み込む方法を示す図である。 現在の時間枠の準最適状態関数を得るための前に計算された状態関数の修正を示す図である。 唯一のレベルによる計算プロセスのメタ制御を示す図である。 複雑な計算プロセスの複数のレベルによるメタ制御および／または階層構造を示す図である。

Claims (35)

1. 最適化システムであって、
   コンピュータ・プログラムを実行するコンピュータ・システムと、
   最適化プログラムと、を有し、
   前記最適化プログラムは、
   最適化問題を受け取り、
   連続型反復パラメータにより前記最適化問題をパラメータ化して連続最適化問題を生成し、
   前記連続最適化問題の解法を採用して前記最適化問題の最適解に向かう軌跡を定め、
   前記最適化問題の複数のパラメータを制御変数として定式化し、前記制御変数を計算するための反復法により最適な軌跡を含む所定領域内で準最適軌跡を求め、
   前記準最適軌跡を繰り返し計算して大域的な最適解を含む所定区間内の準最適解を得る、
   最適化システム。
2. 前記最適化プログラムが、記号形式で表現された最適化問題を受け取る、請求項１に記載の最適化システム。
3. 前記最適化プログラムが、前記記号形式で表現された最適化問題を受け取った後で、
   前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
4. 前記最適化プログラムが、
   前記受け取った最適化問題内の非実数の決定変数を実数の決定変数に変換することによって、前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
5. 前記受け取った最適化問題内の非実数の決定変数を実数の決定変数に変換することが、
   浮動小数点変数を強制的に２つの値０または１のいずれかにするための追加の制約ｙ_ｉ（１−ｙ_ｉ）＝０を伴って範囲［０，１］において、ブール型の各決定変数ｂ_ｉを等価な浮動小数点値ｙ_ｉに変換することと、
   整数型の各決定変数を、等価なｋ項の２進多項式に変換することを含み、ここでｋは前記整数型の決定変数がとり得る最大値の対数（基底は２）にほぼ同等であり、
   さらに、前記２進多項式をｋ個の浮動小数点変数ｙ_０，ｙ_１，．．．，ｙ_ｋ−１によるｋ個の項で表すことを含み、前記変数のすべてが範囲［０，１］に含まれ、ｙ_ｉ（１−ｙ_ｉ）＝０の形のｋ個の制約を伴う、
   請求項４に記載の最適化システム。
6. 前記最適化プログラムが、
   実数型の決定変数を浮動小数点型の決定変数に変換することによって、前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
7. 前記最適化プログラムが、
   等式制約を不等式制約に変換することによって前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
8. 前記等式制約を不等式制約に変換することが、
   個々の等式制約ｈ_ｊ（ｘ）について、
   

   

   のように、ｈ_ｊ（ｘ）を２つの不等式制約に変換し、第３の制約ｓ_ｊ≧０を追加することをさらに含み、ここで、ｍは等式制約を不等式制約に変換する前の不等式制約の数、ｓ_ｊは束縛変数である、請求項７に記載の最適化システム。
9. 前記最適化プログラムが、
   エンベロープ変数を追加すること、および
   前記最適化問題の最適解に向かう軌跡への制約の影響を表すバリア関数を定式化することによって、
   前記受け取った最適化問題から標準形式による前記最適化問題の記号表現を生成する、請求項２に記載の最適化システム。
10. 前記エンベロープ変数を追加することが、
    不等式制約を伴うエンベロープ変数ｚを、等式制約を不等式制約に変換することによって発生する前記最適化問題
    

    

    に加えることをさらに含み、
    前記変数ｚは、
    

    

    として定義され、前記不等式制約は、
    

    

    として定義される、請求項９に記載の最適化システム。
11. 前記バリア関数を追加することが、各制約ｇ_ｉについてバリア関数−ｕ_ｉＢ(ｇ_ｉ)を追加することをさらに含む、請求項９に記載の最適化システム。
12. 前記バリア関数が、各制約ｇ_ｉについて−ｕ_ｉlog(ｇ_ｉ)の形をとる、請求項１１に記載の最適化システム。
13. 前記標準形式による前記最適化問題の記号表現が、
    

    

    と表すことができ、ここで、前記関数Ｆは決定変数ｘ、束縛変数ｓ、エンベロープ変数ｚ、複数の制御変数を備えるバリア関数の係数ベクトルｕ、複数の制御変数を備える束縛変数の係数ベクトルｒの関数であり、ｇ_ｉは不等式制約を定義する式であり、
    前記最適化問題の前記標準形式による記号表現が、さらにコンパクトに、
    

    

    と表すことができ、ここで、
    

    

    であり、ここで、ｍは等式制約を不等式制約に変換する前の不等式制約の数、ｋは変換された不等式制約の数、「状態ベクトル」と呼ばれるｙは決定変数ｘ、束縛変数ｓ、およびエンベロープ変数ｚのベクトル、ベクトルｕはバリア関数の係数ベクトル、ｒは束縛変数の係数ベクトルである、請求項３に記載の最適化システム。
14. 前記最適化プログラムが、
    ｙ_ｋ＋１＝ｙ_ｋ＋φ_ｋ（ｙ_ｋ，ｒ_ｋ，ｕ_ｋ）
    （ｋ＝０，１，..．，ｋ_convergence）
    のように現在の状態ベクトルｙ_ｋから次の状態ベクトルｙ_ｋ＋１を計算できる反復関数φ_ｋをパラメータ化することにより、連続型反復パラメータで前記最適化問題をパラメータ化して連続最適化問題を生成する、請求項１３に記載の最適化システム。
15. 前記最適化プログラムが、０から反復限界Ｔまでの範囲にある連続型反復パラメータτを導入し、ベクトル関数ｙ、ｕ、ｒを連続型反復変数τでパラメータ化することによって、前記反復関数φ_ｋをパラメータ化する、請求項１４に記載の最適化システム。
16. 前記関数Ｆは、前記連続型反復変数τに関してＦ（ｙ（τ），ｕ（τ），ｒ（τ））と表現される、請求項１５に記載の最適化システム。
17. 前記反復関数ｙ_ｋ＋１＝−ｙ_ｋ＋φ_ｋ（ｙ_ｋ，ｒ_ｋ，ｕ_ｋ）が、１次微分方程式
    

    

    に変換され、ここで、前記変数ｙは従属変数であり、ｕおよびｒは制御変数であり、
    前記反復関数がさらに、前記微分方程式の０≦τ≦Ｔにおける解として効率的な最適化の軌跡を計算するための解法を提供する、請求項１６に記載の最適化システム。
18. ニュートン法を使用してφを前記１次微分方程式
    

    

    に変換し、
    ここで、偏導関数ベクトルＧは、
    

    

    として定義され、ヘッシアン行列Ｈは、
    

    

    として定義されるサイズＮｘＮの行列である、請求項１７に記載の最適化システム。
19. 前記最適化プログラムが、前記最適化プログラムのパフォーマンスを最適化するための制御変数を計算する、請求項１７に記載の最適化システム。
20. 前記最適化プログラムが、
    前記制御変数を最適化するための最適化問題を定式化すること、
    前記最適化変数を計算するための反復を定式化すること、および
    第２の連続型反復変数を使用した反復を継続することによって、
    前記最適化プログラムのパフォーマンスを最適化するための制御変数を計算する、請求項１９に記載の最適化システム。
21. 前記最適化プログラムが、連続最適化問題を常微分方程式に変換することによって前記最適化問題の最小の解に向かう軌跡を求める前記連続最適化問題の解法を採用し、
    前記常微分方程式は、連続型反復変数τに関する前記状態ベクトルｙの正規化された変化率を表す以下に示す式であり、
    

    

    ここで、偏導関数ベクトルＧは、
    

    

    として定義され、ヘッシアン行列Ｈは、
    

    

    として定義されるサイズＮｘＮの行列であり、
    前記常微分方程式は、代替的に、
    

    

    と表され、ここで、
    ｗ(τ)＝Ｈ^−１Ｇ(ｙ(τ),ｕ(τ),ｒ(τ))、および
    ρ(ｙ(τ),ｕ(τ),ｒ(τ))＝Ｇ^ＴＨ^−１Ｇ(ｙ(τ),ｕ(τ),ｒ(τ))
    である、請求項１３に記載の最適化システム。
22. 前記最適化プログラムが、前記最適化問題をｕ(τ)，ｒ(τ)，τ∈［０，Ｔ］に関する制御問題
    

    

    で具体化することによって、準最適方式で軌跡を求める制御変数を計算するための反復法を提供する制御問題に前記最適化問題を変換し、
    ここで、Ｑは対称正定値行列であり、前記制御問題は、以下に示す制約条件
    

    

    に従って、前記目的関数Ｆを最適な状態ベクトル変化の軌跡ｙ^＊(τ)内の局所的最小値にするような最適制御ポリシーｕ^＊(τ)およびｒ^＊(τ)を定める、請求項２１に記載の最適化システム。
23. 最適な状態ベクトル変化の最適化軌跡ｙ^＊(τ)を検出する前記最適制御ポリシーｕ^＊(τ)およびｒ^＊(τ)（τ∈［０，Ｔ］）が、ハミルトニアン
    

    

    から獲得され、
    該ハミルトニアンは、ｕおよびｒの関数としてｕ^＊(τ)およびｒ^＊(τ)で最小値をとり、
    前記状態関数ｙ(τ)および共役関数ｐ(τ)の導関数が、すべてのτ∈［０，Ｔ］について、
    

    

    であり、境界条件、
    

    

    を満たすように、連続型共役関数ｐ(τ)が存在する場合に、ｕ^＊(τ)は最適となる、
    請求項２２に記載の最適化システム。
24. 前記最適化プログラムが、前記最適化問題を、準最適方式で前記軌跡を求める制御変数を計算するための反復法を提供する制御問題に変換し、前記制御問題は、
    

    

    として要約され、ここで、ｙ^＊(τ)は準最適状態ベクトル軌跡、ｙ_０は初期状態ベクトル、ｐ(τ)は共役関数、ｕ^＊(τ)およびｒ^＊(τ)は最適制御ポリシー、Ｑは対称正定値行列、τ∈［０，Ｔ］である、請求項２２に記載の最適化システム。
25. 前記最適化プログラムは、準最適軌跡を繰り返し計算して前記最適化問題の準最適解を得る、請求項２３に記載の最適化システム。
26. 最適化問題の当初の準最適解を使用して、修正および展開された連続型準最適解を生成する修正機能をさらに備える、請求項１に記載の最適化システム。
27. 前記修正機能が、
    前記当初の準最適解を時間的に展開することと、
    最適化問題の前記当初の準最適解が計算されたとき以降に発生したイベントを最適に考慮することと、
    特定の時間区間にわたって前記最適化問題の前記当初の準最適解の特定の状態変数を任意選択で固定することと、
    前記最適化問題の前記当初の準最適解を計算するために当初使用された決定変数の値と制約を任意選択で変更することと、
    最適化問題の前記当初の準最適解に関する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することと、
    によって、最適化問題の前記当初の準最適解を使用して、修正および展開された連続型準最適解を生成する、請求項２６に記載の最適化システム。
28. 前記当初の準最適解を時間的に展開することが、
    時刻ｔ_１−ΔＴ＋Ｔからｔ_１−ΔＴ＋Ｔまで長さΔＴの時間区間にわたって前記当初の準最適解を時刻ｔ_１で展開するために、前記当初の準最適解の前記制御関数を使用し、前記当初の準最適解の前記状態関数の偏微分方程式を積分することをさらに含む、請求項２６に記載の最適化システム。
29. 最適化問題の前記当初の準最適解が計算されたとき以降に発生したイベントを最適に考慮することが、
    時刻ｔ_１−ΔＴから時刻ｔ_１までのΔＴの時間区間に発生した前記イベントのすべてを、時刻ｔ_１で発生した単一イベントとして近似し、前記単一イベントを使用して、前記修正および展開された連続型準最適解を計算するときに前記当初の準最適解の初期条件を変更することをさらに含む、
    請求項２６に記載の最適化システム。
30. 特定の時間区間にわたって前記最適化問題の前記当初の準最適解の特定の状態変数を固定することが、
    時刻ｔ_１から時刻ｔ_１＋Ｔまでの特定の時間区間にわたって固定する状態変数を識別することと、
    識別された状態変数のそれぞれについて、前記修正および展開された準最適解のための前記識別された状態変数について計算された値が、前記最適化問題の前期当初の準最適解について計算された値と同じになるように、前記修正および展開された連続型準最適解を計算することと、
    をさらに含む請求項２６に記載の最適化システム。
31. 前記最適化問題の前記当初の準最適解を計算するために当初使用された決定変数の値および制約を変更することが、
    フェンス区間にわたって変数を固定することをさらに含む、請求項３０に記載の最適化システム。
32. 前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
    時間区間Ｔにわたってフェンスのない変数の変動を最小化することをさらに含む、請求項３１に記載の最適化システム。
33. 前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
    前記最適化問題のノミナル解を時刻ｔ_１＋Ｔから時刻ｔ_１＋Ｔ＋ΔＴまで分析的に展開することをさらに含む、請求項３０に記載の最適化システム。
34. 前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
    修正のために特定の区間に関連する変数の優先順位を付けることをさらに含む、請求項３３に記載の最適化システム。
35. 前記最適化問題の前記当初の準最適解に対する前記修正および展開された準最適解の状態関数の変化を最小化し、混乱を最小化することが、
    増分最適化問題の解を展開する間に基準の制約を変更することをさらに含む、請求項３０に記載の最適化システム。
    
    

Images