(57)【要約】【課題】制御対象の動特性推定において、制御対象の入力に含まれる非線形特性を簡便かつ確実な方法により推定する。【解決手段】リレーフィードバックによりリミットサイクルを発生させ、そのときの位相交差周波数におけるプラントゲインとリレー要素出力のデューティー比を測定する。この実験をリレー出力の値を変えて複数回行い、それぞれの実験により得られたデューティー比と基本波におけるゲインの測定値から、プラントの入力に含まれる非線形特性を推定する。【選択図】図１
2005284428.tif 000002
【特許請求の範囲】
【請求項１】
プラントに対して、設定値を設定し、リレー要素を用いて前記リレー要素の出力信号が操作量となるようにリレーフィードバック系を構成し、前記リレー要素は入力の値によって２種類の出力の値を選択して出力するものであり、前記リレー要素の出力の取る値について少なくとも２種類に組み合わせに対して前記設定値の値を同じにして前記リレーフィードバック系にリミットサイクルを発生させ、それらのリミットサイクル発生時における前記リレー要素の出力の信号のデューティー比および前記プラントの出力に係る信号の基本波の振幅に係る量を計測し、それらのリミットサイクル発生時における前記リレー要素出力の基本波の振幅および直流成分の値、および前記プラントの出力に係る信号の基本波の振幅を基に、前記プラントに含まれる非線形性を推定するリミットサイクルによる非線形プラント同定方式。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、フィードバック制御における制御対象の入出力特性を推定する装置に関するものである。
【背景技術】
【０００２】
フィードバック制御において、制御器を設計したり、制御パラメータを設定するのに、プラント（制御対象）の動特性を推定する必要がある。プラントの動特性を推定する方法としては、リレーフィードバックによるリミットサイクル法が実用的な方法として多く用いられている。従来のリミットサイクル法を用いることにより、プラントの位相交差周波数や位相交差周波数におけるゲインを推定することができたり、もう少し一般的な伝達関数を同定することができる。
【０００３】
しかし、一般にプラントは非線形動的システムであり、従来のリミットサイクルテストにおいてはプラントの線形近似モデルしか推定することができず、プラントの非線形特性を推定することができなかった。
【特許文献１】特開２００３−２４１８０１号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
解決しようとする課題は、プラントの非線形特性を簡便かつ確実な方法で推定することである。プラントの非線形性に関しては、一般的にはその構造の自由度が非常に大きいが、実際に制御する際に大きな問題となることが多い入力に含まれる非線形性の推定を目的とする。非線形プラントを図２に示すように入力に含まれる非線形性と線形動的要素の組み合わせにより近似してモデル化を行い、その際の非線形性部分の特性を推定する。
【課題を解決するための手段】
【０００５】
リレーフィードバックによりリミットサイクルを発生させるようにすると、リレー要素の出力の値は離散的な値を取り、連続的に変化するわけではない。したがって、リレー出力を操作量として非線形プラントに入力すると、プラントの入力部にある非線形要素の出力である仮想的な信号も離散的な値を取ることになる。リレー要素の出力が２値信号であれば、非線形要素の出力である仮想的な信号も２値信号となる。この性質を用いて、リレー要素の出力の三つ以上の値に対してプラントの応答を調べることにより、プラントの入力に含まれる非線形性を推定することができる。具体的には次のようにする。
【０００６】
プラントに対してリミットサイクルが発生するようにリレー要素を用いてフィードバック制御系を形成し、リレー要素の出力の取る値について少なくとも２種類に組み合わせに対してリミットサイクルを発生させる。それらのリミットサイクル発生時におけるリレー要素出力の組み合わせとそのときのリレー要素出力信号のデューティー比から、信号の直流成分に着目して、非線形要素の特性に対するある関係式が得られる。また、そのときのリレー要素出力の組み合わせとそのときのリレー要素出力信号のデューティー比およびプラント出力の振幅の組み合わせから、信号の基本波成分に着目して、非線形要素の特性に対するある関係式が得られる。これら２種類の情報を両方用いることによって、少ない実験回数のデータを用いてプラントの非線形要素の特性を精度良く推定する。
【０００７】
信号の直流成分に着目した非線形要素の特性に対する関係式は次のように導出される。リレーフィードバックにより制御量が設定値付近でリミットサイクルを起こしているとき、リレー要素の出力が取る二つの値の組合せを変えても、制御量の時間平均は大体同じ値になる。したがって、図２における非線形要素の出力信号vの時間平均もほぼ同じになる。この関係を、リレー要素出力が取る値の複数の組みに対して立てることにより、非線形関数に対する関係式が得られる。
【０００８】
リミットサイクルにおける信号の基本波成分に着目した非線形要素の特性に対する関係式は次のように導出される。リレーフィードバックにより制御量が設定値付近でリミットサイクルを起こしているとき、リレー要素の出力が取る二つの値の組合せを変えても、リミットサイクルの周期は大きくは変動しないことが多い。非線形要素の出力信号vの基本波成分の大きさは、その信号が取る二つの値とデューティー比によって決まる。また、制御量の基本波成分の大きさは測定可能である。場合によっては、制御量の振幅でもって制御量の基本波の大きさの近似値とすることができる。リミットサイクルの周波数が不変の場合は、非線形要素の出力信号vの基本波成分の大きさと制御量の基本波成分の大きさの比は不変となるので、この関係を用いることにより、非線形関数に対する関係式が得られる。リレー要素の出力が取る二つの値の組合せを変えることによってリミットサイクルの周波数が変化した場合においては、その周波数の変化によるプラントのゲインの変動を補正することもできる。この補正は、リミットサイクルの基本波の周波数付近においては、プラントのゲインの周波数に対する傾きはおおよそ-40dB/dec.となることを利用する。
【０００９】
上記の方法においては、リレー要素の出力が取る二つの値の組合せを変えても、制御量の時間平均がほぼ同じ値になることを仮定しているので、制御量の時間平均の変動が非線形特性の推定誤差につながってしまう。そこで、リレー要素の出力が取る二つの値の組合せを変えても、制御量の時間平均が変動しないように、リレーフィードバックを行う際に、リレー要素の前段に積分を含む要素を追加してもよい。このことにより、非線形特性の推定精度が向上する。このとき、リレー要素の前段に挿入する動的要素の周波数応答を適切に設定することにより、リレーフィードバックシステムにおけるリレー要素の除いた一巡伝達関数の周波数応答のシェーピングを行うことができ、リミットサイクルの周波数を向上させることができる。リミットサイクルの周波数が高くなると、リミットサイクルテストを行うのに要する時間を短縮することができる。そのような動的要素の例として、PID補償器などがある。
【発明の効果】
【００１０】
本発明の非線形プラント同定器を用いることにより、プラントの入力に含まれる非線形特性を少ない実験回数で簡便かつ確実に推定することができる。推定した非線形特性の逆関数を制御器の出力部分に作用させることにより、入力に強い非線形性を持つプラントに対して、その非線形特性を補償することができ、高性能な制御性能を持つ制御系を構築することができる。
【発明を実施するための最良の形態】
【００１１】
本発明を実施するための最良の形態を実施例を通して述べる。
【実施例１】
【００１２】
図１は本発明第１の実施例における非線形性を推定するための同定実験におけるシステム構成図である。リレーフィードバックによるリミットサイクルを発生させる。設定値SVは同定対象であるプラントに対して、どの状態周りでの特性を推定するかを指定する値であり、実際に制御を行う代表的な値であることが望ましい。リレー要素１は入力の信号の符号に応じてuLまたはuHの値を出力するものである。この出力値uLおよびuHは、非線形性推定器３により変化させることができる。リレー出力は操作量u(t)として同定対象であるプラント２に入力される。プラント出力である制御量y(t)はフィードバックされ、設定値SVとの差を演算され、リレー要素１に入力される。リレー要素１の出力の値uLおよびuHを適切に設定すると、リレーフィードバックシステムは自励振動を発生し、図３に示すようなリミットサイクルが生じる。プラントの非線形性を推定するために、リレー要素１の出力の値uLおよびuHを変え、複数回このリレーフィードバックテストを行う。
【００１３】
図２は、本発明第１の実施例で用いる非線形なプラント２に対するモデルである。静的な入力非線形要素２１と線形動的要素２２の組み合わせにより、プラントモデル２０を構成しており、Hammerstein型モデルと呼ばれている。入力非線形要素２１の静特性および線形動的要素２２の動特性を推定する必要がある。線形動的要素２２の動特性に関しては、特開２００３−２４１８０１号公報記載のリミットサイクルによるプラント同定手法を用いて同定を行う。この手法を用いることにより、プラントの入力部における非線形性の影響を受けずに同定することができる。残る問題は入力非線形要素２１の静特性の推定である。プラントモデル２０は一般的な非線形モデルとはなっていないが、実際にはプラントの操作部に含まれる非線形性が顕著であることが多いことなどから、コントローラにおいてプラントの入力部の非線形性を補償するだけで大きな効果が得られることが期待できる。
【００１４】
次に、入力非線形要素２１の特性推定の原理を述べる。この特性推定には、信号の直流成分に注目する方法と、リミットサイクルの基本波成分に注目する方法の二つを併用する。二つの方法を併用することにより、少ない実験回数で入力非線形要素２１の特性推定を精度良く行うことができる。
【００１５】
信号の直流成分に注目する特性推定方法の原理は次のようである。リミットサイクルテストにおいて、制御量の変動の振幅がその定常値に比べて十分に小さい場合や、プラントの位相交差周波数におけるプラントのゲインに対してプラントの直流ゲインが十分に大きい場合においては、入力非線形要素２１の時間平均が設定値SVによりほぼ決まり、リレー要素１の出力の値uLおよびuHにはほとんど依存しない。
【００１６】
いま、リレー要素１の出力の値としてuLkおよびuHk（k=1,・・・）の値をとったときの操作量u(k)の波形が図３に示すようであったとする。そして、そのときのパルス幅pkおよびqkを図のように定義する。また、入力非線形要素２１の入力信号uと出力信号vとの関係をv=f(u)とする。すると、入力非線形要素２１の出力信号vの時間平均が設定値SVの値によりほぼ決まってしまうことになるので、設定値SVが一定のときは、次の関係式が得られる。
【００１７】
（数１）
(pk*f(uHk)+qk*f(uLk))/(pk+qk)=vr（一定）
（k=1,・・・）
【００１８】
ただし、vrは一定の値であり、定常な操作量を入力して定常な制御量の値が設定値SVになるときの操作量に対応する値である。
【００１９】
信号の基本波成分に注目する特性推定方法の原理は次のようである。図１に示すリレーフィードバックにおけるリミットサイクルテストにおいては、リミットサイクルの周波数はプラントの中の線形動的部分であるG(s)の位相交差周波数に近い値となり、リレー要素１の出力の値uLおよびuHにはあまり依存しない。G(s)の位相交差角周波数をωpとする。いま、G(s)の入力である信号v(t)の基本波成分の振幅をak、制御量y(t)の基本波成分の振幅をekとすると、リミットサイクルの角周波数がωpに近いとき、次の近似式が成り立つ。
【００２０】
（数２）
ek ≒ |G(ωp)|*(ωp/ωk)2*ak
【００２１】
（数３）
ωk = 2π/(pk+qk)
【００２２】
各信号の波形の特徴量を図３のように定義すると、信号v(t)の基本波成分の振幅akは次のように求めることができる。
【００２３】
（数４）
ak = 2*(f(uLk)-f(uHk))*sin(π*pk/(pk+qk))/π
【００２４】
また、制御量y(t)の基本波の振幅ekについてであるが、プラントのローパスフィルタ効果により制御量y(t)には高調波成分は多くは含まれておらず、制御量y(t)はほぼ正弦波として扱うことができる。したがって、制御量y(t)の基本波の振幅は、制御量y(t)の振動幅bkを用いて次のように近似することができる。
【００２５】
（数５）
ek = bk/2
【００２６】
実際には、なるべく少ない回数のリミットサイクルテストにより入力非線形要素２１の特性推定を行うことが望ましいので、２回のリミットサイクルテストを行うことにより入力非線形要素２１の特性推定を以下のように行う。
【００２７】
仮定として、操作量u(t)の取る値の範囲は0から1までであるとする。先ず、リレー要素１の出力の取る値をuH1=1、uL1=0と設定し、リミットサイクルテストを行う。実験データから、操作量u(t)のパルス幅p1およびq1、制御量y(t)の振動のピーク・ツー・ピークの値e1を取得する。入力非線形要素２１の入出力特性に対して、f(uL1)=uL1、f(uH1)=uH1とすることにより、ゲインとバイアスに対する冗長性を排除する。そして、次式によりvrの値を決定する。
【００２８】
（数６）
vr = (p1*f(uH1)+q1*f(uL1))/(p1+q1)
【００２９】
さらに、このときのリミットサイクルの各周波数ω1
【００３０】
（数７）
ω1 = 2π/(p1+q1)
【００３１】
における線形要素２２のゲインを求めておく。
【００３２】
（数８）
|G(ω1)| = e1 / a1
【００３３】
（数９）
a1 = 2*(f(uL1)-f(uH1))*sin(π*p1/(p1+q1))/π
【００３４】
（数１０）
e1 = b1/2
【００３５】
次に、リレー要素１の出力の取る値をuH2=0.8、uL2=0.2と設定し、リミットサイクルテストを行う。実験データから、操作量u(t)のパルス幅p2およびq2、制御量y(t)の振動のピーク・ツー・ピークの値e2を取得する。この結果から、次の関係式を得る。
【００３６】
（数１１）
(p2*f(uH2)+q2* f(uL2)) = (p2+q2) * vr
【００３７】
また、
【００３８】
（数１２）
|G(ω2)| = e2 / a2
【００３９】
（数１３）
a2 = 2*(f(uL2)-f(uH2))*sin(π*p2/(p2+q2))/π
【００４０】
（数１４）
e2 = b2/2
【００４１】
の関係および
【００４２】
（数１５）
|G(ω2)| ≒ |G(ω1)|*(ω1/ω2)2
【００４３】
の関係から、次式を得る。
【００４４】
（数１６）
2*(f(uL2)-f(uH2))*sin(π*p2/(p2+q2))/π ≒ e2*(ω2/ω1)2 / |G(ω1)|
【００４５】
数１１および数１６からなる連立方程式を解くことにより、f(uL2)およびf(uH2)の値を求めることができる。
【００４６】
上述の方法により、関数f(u)について、 f(uL1)、f(uH1)、 f(uL2)およびf(uH2)の４点の値が決定できたので、これらの情報をもとに関数f(u)を決定する（図４）。多くの場合、関数f(u)は単調関数であるので、単調関数となるように関数f(u)を決定する。方法としては、関数f(u)の形の候補として、調整パラメータを持ついくつかの関数形を用意しておき、予め決めておいた関数形の順番により、パラメータ調整による関数のあてはめを行い、適切な関数のあてはめができた時点で関数f(u)を決定する。
【００４７】
まず最初に適用する関数系は２次のべき関数である。関数f(u)が上または下に凸である場合に有効な関数形である。f(0)=0、f(1)=1の拘束条件から、関数の形としては次のものを用いる。
【００４８】
（数１７）
f(u) = u + α*u*(1-u)
【００４９】
関数f(u)が0＜u＜1の範囲で単調関数であるための必要十分条件は、-1＜α＜1となる。f(uL2)およびf(uH2)の２点の値がリミットサイクルテストによって得られた値に近くなるようにαの値を決定する。αの値の決定に対しては、最小２乗法を用いてもよいし、直接探索を行ってもよい。得られたαの値に対して、α＜-1となったりα＞1となった場合や、f(uL2)またはf(uH2)の値がリミットサイクルテストによって得られた値と大きく離れてしまった場合においては、２次のべき関数の適用が不適切であると判断して、次のステップに進む。
【００５０】
次に、適用を考える関数形は、次のものである。
【００５１】
（数１８）
f(u) = uα
【００５２】
ここで、αは調整パラメータであり、α＞0である。0＜α＜1に対して関数f(u)は上に凸となり、1＜αに対して下に凸となる。この関数もf(0)=0、f(1)=1となるので、f(uL2)およびf(uH2)の２点の値がリミットサイクルテストによって得られた値に近くなるようにαの値を直接探索により決定する。f(uL2)またはf(uH2)の値がリミットサイクルテストによって得られた値と大きく離れてしまった場合においては、この関数形の適用が不適切であると判断して、次のステップに進む。
【００５３】
次に、３次関数により関数f(u)のあてはめを行うことを検討する。４点のデータを用いて３次関数により関数のあてはめを行うので、３次関数は一意に決まる。ただし、0＜u＜1の範囲において、関数f(u)が単調とならなかった場合においては、３次関数の適用が不適切であると判断して、次のステップに進む。
【００５４】
以上により関数のあてはめが適切にできなかった場合は、折れ線関数を用いる。すなわち、f(u)に対して２点(uL2, f(uL2))、(uH2, f(uH2))を結ぶ直線とし、その値が1を超えていたら1とし、0よりも小さい値であった場合には0とする。
【００５５】
図５は、非線形プラントに対する制御系の一構成である。本発明によりプラント２の非線形要素２１の入出力特性f(u)を推定する。そして、コントローラ４の中に入力非線形補償要素４２を入れ、その入出力特性を推定した関数f(u)の逆関数に設定することにより、入力非線形補償要素４２と非線形要素２１の非線形特性が互いに相殺され、制御演算部から見た制御対象の特性が線形動的要素２２の特性であるG(s)となる。したがって、制御演算部４１は、線形な制御対象に対して設計すればよいことになる。
【００５６】
本発明の第１の実施例においては、同定実験を２種類としてk=1,２としたが、同定実験の種類をさらに多くしてもよい。
【００５７】
本発明の第１の実施例においては、関数f(u)の推定を行っていたが、その逆関数を直接同定するようにしてもよい。コントローラ設計の際に必要となるのは逆関数の方であり、関数f(u)は直接には必要ない。また、関数f(u)の逆関数を推定する際にも、逆関数の単調性を保証するように推定してもよい。
【００５８】
本発明の第１の実施例においては、制御量y(t)の基本波の振幅を制御量y(t)の振動のピーク・ツー・ピークの値を用いて近似していたが、他の方法により制御量y(t)の基本波の振幅を求めてもよい。制御量y(t)に対してsin関数やcos関数との内積を計算して求めてもよい。また、基本波のパワーを算出することにより基本波の振幅を求めてもよい。このとき、制御量y(t)の第４高調波以上の高調波成分は十分に小さいものと仮定し、制御量y(t)の全パワーから第２高調波のパワーを差し引いたものを基本波のパワーとしてもよいし、さらに第３高調波のパワーを差し引いたものとしてもよい。このとき、第２高調波および第３高調波のパワーは、制御量y(t)と矩形波との内積により推定が可能であるので、少ない計算量で制御量y(t)の基本波の振幅を精度良く見積もることができる。
【実施例２】
【００５９】
本発明第２の実施例は、本発明第１の実施例と、プラントモデル、入力非線形要素２１の入出力特性である関数f(u)の推定方法は基本的に同じであるが、リミットサイクルテストにおけるリレーフィードバック系の構成が異なる。
【００６０】
図６は、本発明第２の実施例でのリミットサイクルテストにおけるリレーフィードバック系の構成を示す図である。本発明第１の実施例と比して、リレー要素１の前にＰＩＤコントローラ５が挿入されているところが相異している。ＰＩＤコントローラ５に積分動作を入れることにより、リレー要素１の二つの出力の値を変えても、リミットサイクルテストにおける制御量y(t)の時間平均が全く変化しないので、関数f(u)を精度良く推定することができる。
【００６１】
本発明第２の実施例における関数f(u)の推定方法は、基本的には本発明第１の実施例と同じであるが、信号y(t)の基本波の振幅をＰＩＤコントローラ５の出力であるw(t)に置き換え、かつ、G(s)をG(s)とＰＩＤコントローラ５の伝達関数の積に置き換えるところが相異する。具体的には、次のようである。
【００６２】
実際には、なるべく少ない回数のリミットサイクルテストにより入力非線形要素２１の特性推定を行うことが望ましいので、２回のリミットサイクルテストを行うことにより入力非線形要素２１の特性推定を以下のように行う。
【００６３】
仮定として、操作量u(t)の取る値の範囲は0から1までであるとする。先ず、リレー要素１の出力の取る値をuH1=1、uL1=0と設定し、リミットサイクルテストを行う。実験データから、操作量u(t)のパルス幅p1およびq1、制御量y(t)の振動のピーク・ツー・ピークの値e1を取得する。入力非線形要素２１の入出力特性に対して、f(uL1)=uL1、f(uH1)=uH1とすることにより、ゲインとバイアスに対する冗長性を排除する。そして、次式によりvrの値を決定する。
【００６４】
（数１９）
vr = (p1*f(uH1)+q1*f(uL1))/(p1+q1)
【００６５】
さらに、このときのリミットサイクルの各周波数ω1
【００６６】
（数２０）
ω1 = 2π/(p1+q1)
【００６７】
における線形要素２２とＰＩＤコントローラ５の直列システムのゲインを求めておく。
【００６８】
（数２１）
| C(ω1) G(ω1)| = d1 / a1
【００６９】
（数２２）
a1 = 2*(f(uL1)-f(uH1))*sin(π*p1/(p1+q1))/π
【００７０】
ただし、ＰＩＤコントローラ５の伝達関数をC(s)とした。また、d1はＰＩＤコントローラ５の出力信号v(t)の基本波の振幅である。ここでは、d1の値を、このときの信号v(t)の振動のピーク・ツー・ピークの値c1を用いて次のように評価する。
【００７１】
（数２３）
d1 = c1/2
【００７２】
次に、リレー要素１の出力の取る値をuH2=0.8、uL2=0.2と設定し、リミットサイクルテストを行う。実験データから、操作量u(t)のパルス幅p2およびq2、信号v(t)の振動のピーク・ツー・ピークの値c2を取得する。この結果から、次の関係式を得る。
【００７３】
（数２４）
(p2*f(uH2)+q2* f(uL2)) = (p2+q2) * vr
【００７４】
また、
【００７５】
（数２５）
|C(ω2) G(ω2)| = d2 / a2
【００７６】
（数２６）
a2 = 2*(f(uL2)-f(uH2))*sin(π*p2/(p2+q2))/π
【００７７】
（数２７）
d2 = c2/2
【００７８】
の関係および
【００７９】
（数２８）
| C(ω2)G(ω2)| ≒ | C(ω1)G(ω1)|*(ω1/ω2)2
【００８０】
の関係から、次式を得る。
【００８１】
（数２９）
2*(f(uL2)-f(uH2))*sin(π*p2/(p2+q2))/π ≒ e2*(ω2/ω1)2 / | C(ω1)G(ω1)|
【００８２】
数２４および数２９からなる連立方程式を解くことにより、f(uL2)およびf(uH2)の値を求めることができる。
【００８３】
上述の方法により、関数f(u)について、 f(uL1)、f(uH1)、 f(uL2)およびf(uH2)の４点の値が決定できたので、これらの情報をもとに関数f(u)を決定する（図４）。この方法は本発明第１の実施例と同じである。
【００８４】
本発明第２の実施例における本発明第２の実施例に対する他の優位な点として、ＰＩＤコントローラ５のパラメータを適切に設定することにより、リミットサイクルの周期を短くすることができ、リミットサイクルテストに要する時間を短縮することが可能であることが挙げられる。ＰＩＤコントローラ５の微分動作により、線形要素２２とＰＩＤコントローラ５を直列にしたシステムの周波数シェーピングが可能となるためである。
【産業上の利用可能性】
【００８５】
本発明のリミットサイクルによる非線形プラント同定器は、制御系設計の前段階に行うプラント同定において利用することができる。また、プラントの特性推定機能を持つコントローラに実装することにより、非線形プラントに対応したオートチューニング機能つき非線形コントローラを実現することができる。
【図面の簡単な説明】
【００８６】
【図１】本発明第１の実施例における同定実験におけるシステム構成図。
【図２】本発明において用いる非線形プラントのモデル構造。
【図３】リレーフィードバックによるリミットサイクルにおける操作量と制御量の信号波形の説明図。
【図４】プラントの入力部における非線形特性の説明の図。
【図５】非線形特性補償機能付きコントローラおよび制御系の説明図。
【図６】本発明第２の実施例における同定実験におけるシステム構成図。
【符号の説明】
【００８７】
１・・・リレー要素
２・・・プラント
２０・・・プラントモデル
２１・・・入力非線形要素
２２・・・線形動的要素
３・・・非線形性推定器
４・・・コントローラ
４１・・・制御演算部
４２・・・入力非線形性補償要素
５・・・ＰＩＤコントローラ
【図１】
2005284428.tif 000003
【図２】
2005284428.tif 000004
【図３】
2005284428.tif 000005
【図４】
2005284428.tif 000006
【図５】
2005284428.tif 000007
【図６】
2005284428.tif 000008