(57)【要約】本発明は技術システムを設計する方法に関する。本方法では、技術システムはパラメータに依存する所定数の目的関数によってモデル化される。その際、個々の目的関数は重み係数で重み付けされる。本方法は前記パラメータと前記重み係数を変数として有する方程式系を変数空間内で解く。なお、この方程式系の解は変数空間内で解空間の作動点を形成している。本発明による方法では、作動点は予測子修正子法により決定される。その際、第１の作動点に基づいて、確率量から形成される予測子が変数空間内で求められ、続いて修正ステップにおいて第２の作動点が求められる。求められた作動点は技術システムの設計に使用される。技術システムの設計は既存の技術システムの新たな設計又は変更ないしは適合化であってもよい。
2005537539.tif 000002
【特許請求の範囲】
【請求項１】
技術システムを設計する方法において、
技術システムはパラメータに依存する所定数の目的関数を有し、
前記目的関数の各々を重み係数で重み付けし、
前記パラメータと前記重み係数を変数として有する方程式系を変数空間内で解き、
前記方程式系の解は前記変数空間内で解空間の作動点を形成し、
前記作動点を予測子修正子法により決定し、その際、確率量によって形成される予測子を第１の作動点に基づいて前記変数空間内で決定し、続いて修正ステップにおいて第２の作動点を求め、
前記作動点を技術システムの設計に使用することを特徴とする、技術システムを設計する方法。
【請求項２】
前記予測子を乱数によって決定する、請求項１記載の方法。
【請求項３】
前記確率変数は正規分布する、請求項１又は２記載の方法。
【請求項４】
前記確率量は下記の方程式を満たす確率過程Ｚｔであり、
ｄＺｔ＝εＰ（Ｚｔ）ｄＢｔ
ここで、
Ｐ（Ｚ）は有効な作動点ｚにおける解空間の接空間への射影行列であり、
εはスケール係数であり、
【数１】
2005537539.tif 000008
は変数空間内でのブラウン運動である、請求項１から３のいずれか１項記載の方法。
【請求項５】
前記作動点としてパレート最適な点を求める、請求項１から４のいずれか１項記載の方法。
【請求項６】
前記作動点は解空間内の正の重み係数を有する点である、請求項１から５のいずれか１項記載の方法。
【請求項７】
前記作動点は１つ又は複数の制約条件を満たし、ただし、前記１つ又は複数の制約条件は方程式系の別の変数によって変数空間内に表現される、請求項１から６のいずれか１項記載の方法。
【請求項８】
前記制約条件は等号制約条件及び／又は不等号制約条件である、請求項７記載の方法。
【請求項９】
前記不等号制約条件をスラック変数を用いて方程式制約条件に変換する、請求項８記載の方法。
【請求項１０】
前記解空間は多様体、とりわけ変数空間内の部分多様体である、請求項１から９のいずれか１項記載の方法。
【請求項１１】
前記第１の作動点を重み付け法により決定する、請求項１から１０のいずれか１項記載の方法。
【請求項１２】
前記第１の作動点において前記解空間に対する接平面を求め、該接平面内で前記予測子を確定する、請求項１から１１のいずれか１項記載の方法。
【請求項１３】
１つ又は複数の負の重み係数を有する負の予測子が現れた場合、作動点の解空間の部分平面での鏡映によって新たな予測子を求める、請求項１から１２のいずれか１項記載の方法。
【請求項１４】
前記第１の作動点と前記負の予測子との間を走る軌跡と前記解空間の部分平面との交点を求め、
前記交点と前記負の予測子とにより張られたベクトルの、前記解空間の部分平面に対する接線成分を求め、ただし、前記部分平面の点においては、負の予測子に対しては負であった前記重み係数はゼロに等しく、
前記交点と前記負の予測子とにより張られたベクトルの前記接線成分に対する法線成分を求め、
前記法線成分を前記負の予測子から２回減算することにより新たな予測子を求める、請求項１３記載の方法。
【請求項１５】
前記修正子法はニュートン法である、請求項１から１４のいずれか１項記載の方法。
【請求項１６】
前記作動点を予測子修正子法の反復により求め、その際、反復ステップにおいて、予測子修正子法の第１の作動点として、先行する反復ステップの第２の作動点を使用する、請求項１から１５のいずれか１項記載の方法。
【請求項１７】
前記方法を中断条件により終了する、請求項１６記載の方法。
【請求項１８】
前記中断条件は、所定数の作動点が求められ、及び／又は所定の制限時間に達したときに満たされる、請求項１７記載の方法。
【請求項１９】
請求項１から１８のいずれか１項記載の方法を実行する技術システムを設計するための装置において、
該装置はプロセッサユニットを有しており、
該プロセッサユニットは、
前記技術システムがパラメータに依存する所定数の目的関数を有し、
前記目的関数の各々が重み係数で重み付けされ、
前記パラメータと前記重み係数を変数として有する方程式系が変数空間内で可解であり、
前記方程式系の解が前記変数空間内で解空間の作動点を形成し、
前記作動点が予測子修正子法により決定可能であり、ただし、確率量によって形成される予測子が第１の作動点に基づいて前記変数空間内で決定され、続いて修正ステップにおいて第２の作動点が求められ、
前記作動点が技術システムの設計に使用可能であるように構成されている、ことを特徴とする技術システムを設計するための装置。
【請求項２０】
確立量を生成するための乱数発生器を有する、請求項１９記載の装置。
【請求項２１】
計算機上で実行可能であり、請求項１から１８のいずれか１項記載の方法を実行することのできるコンピュータプログラムが記憶された記憶媒体を有することを特徴とするコンピュータプログラム製品。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は技術システムを設計するための方法及び装置に関する。
【０００２】
複雑な技術システムの設計のためには、しばしば複数の互いに対立する基準に関してこのシステムを最適化することが必要となる。これらの基準は例えば製造コストや効率性などのようなシステムの目的関数に影響を与える。さらに、システムの可能な作動点は制約条件によって制限されることがありうる。このことは、システムの最適作動点の集合を求めるという問題、すなわち、これらの作動点をすべての基準に関してそれ以上同時に最適化することが不可能であるようなシステムの可能な作動点の集合を決定するという問題を生じる。この最適点の集合から、ユーザは隠された基準又は専門知識を考慮して使用に最も適したシステム作動点を選択することができる。
【０００３】
［１］から、複数の基準に関して技術システムを最適化する重み付け法が公知である。この方法では、スケールパラメータによるスカラー値最適化問題への変換が使用されている。この方法は、非常に多くのスカラー値最適化を実行しなければならないため、数値計算的に非常にコスト高であるという欠点を有している。さらに、スケールパラメータの選択及び変更はユーザとのインタラクションを必要とし、これに関しては自動化ができない。
【０００４】
［２］には、複数の基準に関して技術システムを最適化する確率的方法が記載されている。この方法では、最適化問題を解くために確率微分方程式が使用される。この方法は、多数の２乗最適化問題を解かなくてはならないため、数値計算的に非常にコスト高であるという欠点を有している。さらに別の欠点は、この方法では個々の目的関数が重み付けされず、そのためユーザが最適点の選択のために重要な情報を使用することができないことにある。
【０００５】
［３］からは、複数の基準に関して技術システムを最適化するホモトピー法が公知である。この方法では、目的関数に対する重み係数の他に、制約条件を考慮するためにラグランジュ乗数が使用される。この方法の欠点は、ユーザとのインタラクションが必要であり、これに関してこの方法を自動化することができないということにある。
【０００６】
したがって、本発明の課題は数値計算的に効率的な自動化された技術システム設計方法を提供することである。
【０００７】
この課題は独立請求項に記載さている特徴によって解決される。本発明の発展形態は従属請求項から得られる。
【０００８】
本発明による技術システムを設計する方法では、技術システムはパラメータに依存する所定数の目的関数によってモデル化される。その際、個々の目的関数は重み係数で重み付けされる。本方法は前記パラメータと前記重み係数を変数として有する方程式系を変数空間内で解く。なお、この方程式系の解は変数空間内で解空間の作動点を形成する。本発明による方法では、作動点は予測子修正子法により決定される。その際、確率量から形成される予測子が第１の作動点に基づいて変数空間内で求められ、続いて修正ステップにおいて第２の作動点が求められる。求められた作動点は技術システムの設計に使用される。技術システムの設計は既存の技術システムの新たな設計又は変更ないしは適合化であってもよい。
【０００９】
本発明の利点は、確率量による予測子の形成によってこの方法が自動化され、ユーザの介入をもはや必要としないということにある。数値計算的な予測子修正子法を確率的方法と結びつけることによって、技術システムの作動点の計算のために計算機資源が効率的に利用されることが保証される。
【００１０】
本発明の１つの有利な実施形態では、予測子は乱数によって定められるので、とりわけ方法実行の際に乱数発生器を使用することができ、そのため容易に方法の自動化が保証される。
【００１１】
本発明の特に有利な別の実施形態では、乱数は正規分布する。これにより、方法の進行に伴って解空間内に形成される作動点の軌跡が解空間全体において一様分布を、したがって可能なすべての作動点の良好な被覆を保証する。正規分布する乱数を使用することにより、本発明による方法によって特に解空間上のブラウン運動をモデル化することができる。
【００１２】
有利には、本発明によって求められる作動点は、すべての目的関数に関してそれ以上最適化が不可能ないわゆるパレート最適点である。とりわけ、本発明による方法では、作動点として、正の重み係数を有する解空間内の点が選択される。
【００１３】
本発明のさらに別の有利な実施形態では、作動点はさらに１つ又は複数の制約条件を満たさなければならない。なお、各々の制約条件は方程式系の別の変数によって変数空間内に表現される。これらの制約条件は等号制約条件及び／又は不等号制約条件であってよい。不等号制約条件の場合には、有利にはスラック変数が導入され、このスラック変数によって不等号制約条件を等号制約条件に変換することができる。スラック変数の利用については実施例の詳細な説明においてより詳しく説明する。
【００１４】
有利には、作動点の解空間は多様体、とりわけ変数空間内の部分多様体である。［３］においては、どのような前提条件の下で解空間がそのような多様体を形成するかが説明されている。
【００１５】
とりわけ本発明による方法の開始時には差し当たり第１の作動点が存在しているので、本発明の特別な実施形態では、この第１の有効な作動点が重み付け法により決定される。なお、重み付け法の使用はすでに従来技術から公知である（［１］参照）。
【００１６】
本発明で使用される予測子修正子法では、有利には、第１の作動点において解空間への接平面が求められ、続いてこの接平面内で予測子が確定される。
【００１７】
本発明による方法の１つの発展形態においては、１つ又は複数の負の重み係数を有する負の予測子が生じた場合には、有効な作動点の解空間の部分平面での鏡映によって新たな予測子が決定される。これにより、有効な作動点の新たな領域を求めることができる。なお、これらの作動点は、隠された付加的基準又はユーザの専門知識との関連でユーザにとって特別に重要であることがありうる。
【００１８】
本発明の１つの有利な実施形態では、鏡映ステップにおいて、第１の作動点と負の予測子との間を走る軌跡と解空間の部分平面との交点が求められる。そして、交点と負の予測子とにより張られたベクトルの、解空間の当該部分平面への接線成分が求められる。ただし、負の予測子に対しては負であった重み係数が今やこの部分平面の点ではゼロに等しい。続いて、交点と負の予測子とにより張られたベクトルのこの接線成分に対する法線成分が求められる。最後に、法線成分を負の予測子から２回減算することにより新たな予測子が求められる。
【００１９】
本発明で使用される修正子法としては、有利には、従来技術から周知のニュートン法が使用される。ニュートン法は数値計算的に容易に実施可能である。
【００２０】
有利には、作動点は予測子修正子法の反復によって決定される。その際、反復ステップにおいて、予測子修正子法の第１の作動点として、先行する反復ステップの第２の作動点が使用される。この方法は例えば中断条件によって終了される。中断条件は、本発明の有利な実施形態においては、所定数の作動点が求められ、及び／又は所定の制限時間に達したときに満たされる。
【００２１】
上述した技術システムを設計する方法の他に、本発明は上記方法を実行することのできる技術システムを設計するための装置にも関している。とりわけ、この方法は確率量によって予測子を形成することのできるプロセッサユニットを有している。
【００２２】
有利には、この装置は確率量を表す乱数を生成するための乱数発生器を有している。
【００２３】
また、本発明又は上述した各々の発展形態は、本発明を実施する計算機により実行可能なコンピュータプログラムが記憶された記憶媒体を有するコンピュータプログラム製品によって実現することもできる。
【００２４】
以下では、図に基づいて本発明の実施例を示し、説明する。
【００２５】
図１は、本発明による技術システム設計方法のフローチャートであり、
図２は、本発明で使用される予測子修正子法を具体的に示したスケッチであり、
図３は、本発明の変化形で使用される鏡映法を具体的に示したスケッチであり、
図４は、本発明による方法を実行するためのプロセッサユニットを示している。
【００２６】
図１には、本発明による技術システム設計方法のフローチャートが示されている。
【００２７】
まず、ステップ１０１において、技術システムの記述形式が選択される。この記述形式は所定数ｋの目的関数ｆ＝（ｆ１，…，ｆｋ）を含んでいる。ただし、これらの目的関数はそれぞれ技術システムのｎ個のプリセットパラメータｘ１〜ｘｎに依存している。これらの目的関数は例えば資本コストｆ１及び効率性ｆ２である。この場合、目的関数は以下の式で表される。
【００２８】
【数１】
2005537539.tif 000003
ここで、ｘ＝（ｘ１，…，ｘｎ）である。
パラメータｘ１〜ｘｎは技術システムの設計パラメータ又は動作パラメータであってよい。
【００２９】
本発明による方法によれば、技術システムの設計に使用される有効な作動点は目的関数をパラメータに関して最適化することにより決定される。ただし、最適化基準は一般に互いに競合するので、必ずしもすべての目的関数ｆ１〜ｆｋが同時に最適化されるわけではない。
【００３０】
技術システムは有効な作動点においてさらに所定数ｍの制約条件ｈ＝（ｈ１（ｘ），…，ｈｍ（ｘ））によって制限される。このことは次の式で表される。
【００３１】
ｈ（ｘ）＝０
ここで、０＝（０，…，０）はゼロベクトルを表している。上式は等号制約条件であり、択一的に不等号制約条件を考慮してもよい。このような不等号制約条件は例えば次のようである。
【００３２】
ｈ（ｘ）＜０  又は  ｈ（ｘ）＞０
不等号制約条件付き最適化問題を解くために、ｍ個のスラック変数ｓ＝（ｓ１，…，ｓｍ）が導入される。不等号制約条件はスラック変数ｓにより以下の等号制約条件に変形することができる。
【００３３】
ｈ（ｘ）＋ｓ＝０  又は  ｈ（ｘ）−ｓ＝０
本発明のこの実施形態において説明される最適化法では、有効な作動点とは以下の条件を満たすいわゆるパレート最適な点である。
【００３４】
【数２】
2005537539.tif 000004
この最適化問題の解は次の非線形方程式系の解であることが示すことができる。
【００３５】
【数３】
2005537539.tif 000005
ここでは、制約条件がラグランジュ乗数λ＝（λ１，…，λｍ）によって考慮され、目的関数ｆｉは重み係数αｉで重み付けされる。ただし、すべての重み係数の和は１に正規化されている、すなわち、
【００３６】
【数４】
2005537539.tif 000006
であることに注意しなければならない。とりわけ、個々の重み係数αｉは負又はゼロであってもよい。したがって、最適化問題の解は上記方程式系の（ｎ＋ｍ＋ｋ）次元変数空間内のベクトル（ｘ，λ，α）である。
【００３７】
［３］で示されているように、この方程式系の解は所定の前提条件の下で変数空間内の（ｋ−１）次元部分多様体Ｍを記述している。
【００３８】
以下で説明する有効作動点を求めるための数値計算ステップは、実質的に［３］に記載されているホモトピー法に基づいている。ホモトピー法では、パレート最適点を求めるために予測子修正子法が使用される。
【００３９】
ステップ１０２では、技術システムの記述形式１０１に基づいて、第１のパレート最適点が例えば重み付け法のような標準的方法によって求められる。
【００４０】
次のステップ１０３では、この第１のパレート最適点において、有効作動点の多様体Ｍへの点ｚにおける（ｋ−１）次元接平面ＴｚＭが求められる。このために、点ｚにおいて方程式系Ｆのヤコビ行列にＱＲ分解が施される。このＱＲ分解から、接平面を張る直交基底｛ｑ１…ｑｋ−１｝が決定される。この際に実行される個々の数値計算ステップは［３］に詳細に記載されている。
【００４１】
次のステップ１０４では、この接平面内で予測子ｙが決定される。ただし、予測子は−［３］に記載されているホモトピー法とは異なって−正規分布するｋ−１次元の乱数ベクトルｂによって接平面内に形成される。予測子ｙは以下の形式を有している。
【００４２】
ｙ＝ｚ＋（ｑ１…ｑｋ−１）ｂ
このような乱数ベクトルの使用により、部分多様体Ｍ上のブラウン運動をモデル化することができる。その際、ブラウン運動は近似的に以下のように表される。
【００４３】
ｄＺｔ＝εＰ（Ｚｔ）ｄＢｔ
ここで、
Ｐ（Ｚ）は有効な作動点ｚにおける解空間の接空間への射影行列であり、
εはスケール係数であり、
【００４４】
【数５】
2005537539.tif 000007
は変数空間内でのブラウン運動である。
【００４５】
このブラウン運動のモデル化のために、ｂに対してｋ−１次元の正規分布Ｎ（０ｋ−１，ｔΔεＩｋ−１）が選択される。ここで、平均値は（ｋ−１）次元のゼロベクトル０ｋ−１であり、分散は（ｋ−１）次元の単位行列Ｉｋ−１にブラウン運動の歩幅ｔΔ及びスケール係数εを乗じたものである。
【００４６】
択一的に、まず（ｍ＋ｎ＋ｋ）次元変数空間内で正規分布する乱数ベクトルを求め、続いてこの乱数ベクトルを（ｋ−１）次元接平面ＴｚＭに射影することによって、予測子を求めてもよい。
【００４７】
その後で、ステップ１０５において、予測子は修正子法によってパレート最適点の多様体上に射影される。なおここで、修正子法は例えば数値計算におけるニュートン法である。このようにして、新たな有効作動点がパレート最適点の多様体上に求められる。
【００４８】
ステップ１０３，１０４及び１０５は反復実行される。その際、新たな有効作動点の計算のための出発点として、先行する反復ステップにおいて求められた作動点が使用される。
【００４９】
ステップ１０６では、中断条件が満たされているか否か、すなわち、例えば所定数の反復が実行されたか又は所定の制限時間に達したかが調べられる。中断条件が満たされていない場合には、ステップ１０３に戻り、次の反復が実行される。これは中断条件が満たされるまで続けられる。
【００５０】
ステップ１０６において中断条件が満たされている場合には、次のステップ１０７において、特定のパレート最適点の集合が、重み係数αｉが正であるような点に限定される。
【００５１】
最後のステップ１０８では、これらのパレート最適点からユーザの要求にとって効率的な技術システムの作動点がユーザによって選択され、技術システムはこの効率的な作動点を用いて設計される。
【００５２】
図２は、本発明で使用される予測子修正子法の２次元グラフ表示が示されている。図２では、ｚｉが部分多様体Ｍ上のパレート最適点を表している。なお、この点はｉ番目の反復ステップにおいて得られたものである。新たなパレート最適点を求めるために、まず部分多様体Ｍに対する点ｚｉにおける接平面ＴｚｉＭが求められる。この接平面は図２では破線で示されている。次のステップ１０４では、正規分布する乱数を用いて、接平面ＴｚｉＭ内で予測子点ｙｉ＋１が求められる。その次の修正ステップ１０５では、新たなパレート最適点ｚｉ＋１が求められる。なお、修正ステップ１０５は例えばニュートン法であってよい。この方法は続行され、その際、新たな予測ステップの出発点としてパレート最適点ｚｉ＋１が使用される。
【００５３】
図３は、本発明による方法の変化形に関するものである。ここでは、負の重み係数αｉを有する予測子が生じると、正のαｉを有する新たな予測子を求めるために鏡映が行われる。図３はこの射影ステップの実行を３次元表示で示したものである。
【００５４】
図３には、パレート最適点ｚに基づいて、負のαｉを有する予測子ｙｎｅｇが求められるケースが示されている。これは、点ｚと点ｙの間のコースが点Ｓで接平面ＴｚＭを貫通することによってグラフで具体的に示される。点Ｓは接平面ＴｚＭの部分平面上にあり、点Ｓに対する座標αｉは値ゼロを有している。鏡映を実行するためには、まず交点Ｓが求められる。これは、点ｚとｙとを通る直線のパラメータ表示からαｉ成分を射影によって取り出す射影作用素を使用することで行うことができる。点Ｓを求めた後には、点Ｓとｙとの間を走るベクトルｘｎｅｇを求めることができる。そして、このベクトルは部分平面における接線成分ｔと法線成分ｎとに分解される。よって、接線成分に関してｔ＝ｘｎｅｇ−ｎが成り立つ。続いて鏡映ステップが実行される。その際、鏡映された新たなベクトルｘｎｅｕは古いベクトルｘｎｅｇと同じ接線成分ｔを有し、法線成分は符号を逆にした古いベクトルｘｎｅｇの法線成分ｎに相当する。したがって、新たなベクトルに関しては、ｘｎｅｕ＝ｔ−ｎ＝（ｘｎｅｇ−ｎ）−ｎ＝ｘｎｅｇ−２ｎが成り立つ。これにより、接平面ＴｚＭで鏡映された新たな予測子ｙｎｅｕが得られる。上述した鏡映法は方法の数値計算的効率性を高める。というのも、負の重み係数αｉを有する点の形成が防止され、したがって計算技術上の資源がより効率的に利用されるからである。
【００５５】
図４には、本発明による方法を実行するためのプロセッサユニットＰＲＺＥが示されている。プロセッサユニットＰＲＺＥは、プロセッサＣＰＵ、メモリＭＥＭ、及び入出力インタフェースＩＯＳを有しており、入出力インタフェースＩＯＳはインタフェースＩＦＣを介して様々に利用される：グラフィックインタフェースを介して、モニタＭＯＮ上の出力が視覚化される及び／又はプリンタＰＲＴに出力される。入力はマウスＭＡＳ又はキーボードＴＡＳＴを介して行われる。また、プロセッサユニットＰＲＺＥは、メモリＭＥＭとプロセッサＣＰＵと入出力インタフェースＩＯＳとの間の接続を可能にするデータバスＢＵＳを使用することができる。さらに、データバスＢＵＳには、付加的なメモリ、データストレージ（ハードディスク）又はスキャナなどの付加的なコンポーネントを接続することができる。
【図面の簡単な説明】
【００５６】
【図１】本発明による技術システム設計方法のフローチャートを示す。
【図２】本発明で使用される予測子修正子法を具体的に示したスケッチを示す。
【図３】本発明の変化形で使用される鏡映法を具体的に示したスケッチを示す。
【図４】本発明による方法を実行するためのプロセッサユニットを示す。文献一覧表：［１］Ｃ．Ｈｉｌｌｅｒｍｅｉｅｒ：”Ｎｏｎｌｉｎｅａｒ  Ｍｕｌｔｉｏｂｊｅｃｔｉｖｅ  Ｏｐｔｉｍｉｚａｔｉｏｎ：Ａ  Ｇｅｎｅｒａｌｉｚｅｄ  Ａｐｐｒｏａｃｈ”，第３．２章、Ｂｉｒｋｈａｅｕｓｅｒ  Ｖｅｒｌａｇ，２００１［２］Ｃ．Ｈｉｌｌｅｒｍｅｉｅｒ：”Ｎｏｎｌｉｎｅａｒ  Ｍｕｌｔｉｏｂｊｅｃｔｉｖｅ  Ｏｐｔｉｍｉｚａｔｉｏｎ：Ａ  Ｇｅｎｅｒａｌｉｚｅｄ  Ａｐｐｒｏａｃｈ”，第３．３章、Ｂｉｒｋｈａｅｕｓｅｒ  Ｖｅｒｌａｇ，２００１［３］Ｃ．Ｈｉｌｌｅｒｍｅｉｅｒ：”Ｎｏｎｌｉｎｅａｒ  Ｍｕｌｔｉｏｂｊｅｃｔｉｖｅ  Ｏｐｔｉｍｉｚａｔｉｏｎ”，Ｊｏｕｒｎａｌ  ｏｆ  Ｏｐｔｉｍｉｚａｔｉｏｎ  Ｔｈｅｏｒｙ  ａｎｄ  Ａｐｐｌｉｃａｔｉｏｎ，Ｖｏｌ．１１０／３，ｐｐ．５５７−５８３，Ｐｌｅｎｕｍ  Ｐｒｅｓｓ，Ｎｅｗ  Ｙｏｒｋ，２００１
【図１】
2005537539.tif 000009
【図２】
2005537539.tif 000010
【図３】
2005537539.tif 000011
【図４】
2005537539.tif 000012
【国際調査報告】
2005537539.tif 000013
2005537539.tif 000014
2005537539.tif 000015
2005537539.tif 000016
2005537539.tif 000017
2005537539.tif 000018