"設定方位ψRに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、船体パラメータのパラメータ不確かさΔを加味したときの特性多項式DcΔを、以下で表したときに、パラメータ不確かさΔが第1仕様値Δspecのときに、DcΔsの特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζeとして、そのζeが第2仕様値ζspecの値を満足するωe/ωf=ρを求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする船舶用自動操舵装置。前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Dc4Δを、としたときに、Dc4Δsの特性根のうちの最も小さい減衰係数ζをζe2とし、ζe2が第2仕様値ζspecの値を満足するωe2/ωf=ρ2を求め、このρ2を用いた初期値とした数値解法によって前記ωe/ωf=ρを求めることを特徴とする請求項1記載の船舶用自動操舵装置。前記初期値は、で表されることを特徴とする請求項2記載の船舶用自動操舵装置。前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインKsnと時定数Tsnの比の不確かさであることを特徴とする請求項1ないし3のいずれか1項に記載の船舶用自動操舵装置。前記特性多項式DcΔは、で表され、ここで、Dwsは、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、と表され、Bsは、で表されることを特徴とする請求項4記載の船舶用自動操舵装置。しかしながら、従来においては、パラメータ不確かさおよび外乱モデルの影響を間接的、近似的に扱うために、閉ループ安定性に多少の変動を持つという問題がある。即ち、推定ループにおいて閉ループ安定性を劣化させる場合があるが、従来の設計による装置では、推定ループにおけるパラメータ不確かさの影響を直接扱っていない。また、船体モデルはパラメータ不確かさを複数もつために、特許文献2で採用したパラメータ不確かさ以外が船種、海象、載貨などにより無視できない状況を生じる場合があり、用いたパラメータ不確かさによる影響が減衰係数に間接的に作用するために局所的な安定性しか確保することができない場合がある、という問題がある。本発明はかかる課題に鑑みなされたもので、外乱モデルを含んだ推定器と操舵系の閉ループにおいて、パラメータ不確かさに起因する推定誤差を低減することができる推定ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。また、本発明の更なる目的は、推定ループの根を不安定にさせるパラメータ不確かさを特定し、それを用いて制御ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。前述した目的を達成するために、請求項1記載の発明は、設定方位ψRに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、船体パラメータのパラメータ不確かさΔを加味したときの特性多項式DcΔを、以下で表したときに、パラメータ不確かさΔが第1仕様値Δspecのときに、DcΔsの特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζeとして、そのζeが第2仕様値ζspecの値を満足するωe/ωf=ρを求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする。請求項2記載の発明は、請求項1記載の船舶用自動操舵装置において、前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Dc4Δを、としたときに、Dc4Δsの特性根のうちの最も小さい減衰係数ζをζe2とし、ζe2が第2仕様値ζspecの値を満足するωe2/ωf=ρ2を求め、このρ2を用いた初期値とした数値解法によって前記ωe/ωf=ρを求めることを特徴とする。請求項3記載の発明は、請求項2記載の船舶用自動操舵装置において、前記初期値は、で表されることを特徴とする。請求項4記載の発明は、請求項1ないし3のいずれか1項に記載の船舶用自動操舵装置において、前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインKsnと時定数Tsnの比の不確かさであることを特徴とする。請求項5記載の発明は、請求項4記載の船舶用自動操舵装置において、前記特性多項式DcΔは、で表され、ここで、Dwsは、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、と表され、Bsは、で表されることを特徴とする。本発明によれば、外乱モデルとしての波浪モデル及び舵角オフセットモデルを含んだ推定器を含む閉ループの特性多項式に基づき、特性多項式に船体パラメータのパラメータ不確かさΔを加味したときの特性多項式DcΔからパラメータ不確かさΔが第1仕様値Δspecのときに、DcΔsの特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζeとして、そのζeが第2仕様値ζspecの値を満足するωe/ωf=ρを求め、推定器の船体モデルの固有角周波数ωeを設定することから、閉ループ安定性を増加させることができる。また、前記パラメータ不確かさΔを、船体パラメータの旋回力ゲインKsnと時定数Tsnの比の不確かさとすることで、閉ループ安定性に及ぼす影響の高いパラメータ不確かさに起因する推定根の不安定性を排除し、従来のものより信頼性を向上させることができる。外乱モデルを含まない推定器の推定根から出発して、外乱モデルを含んだ推定器の推定根の予測値を求めて、特性多項式DcΔの特性根を求めることで、求解を高い収束性で短時間に行うことができるようになる。以下、図面を用いて本発明の実施の形態を説明する。1.1保針系の構成本発明で対象とする船舶用自動操舵装置10の保針モードは、図1に示すように操舵機14及び船体16を含む制御対象18を設定方位ψRに船首方位ψを追従させるため、フィードバック制御器12で偏差から演算された舵角指令δcによって制御する。δcは操舵機14の舵角δを経て船体に作用する。船体はパラメータ不確かさをもち外乱が印加する。操舵機14は舵角および舵速度の振幅制限をもつが、保針時のδcおよびの振幅が制限以下と見なせるので操舵機を省略するδ=δc。なお設定方位ψRは保針時一定になるため、簡単化のためψR=0とする。1.2制御対象とパラメータ不確かさ図1に示した制御対象18は、操舵機14を除き、船体モデルと外乱モデルとから構成される。1.2.1船体モデル船体モデルは、図1Aに示される野本の線形モデルを用いる。ここで、sはラプラス演算子を、ψは船首方位を、rを旋回角速度を、Ksn,Tsn,は操縦性指数船体パラメータでそれぞれ旋回力ゲイン、時定数のノミナル値添え字s:船体、n:ノミナル値を示す。尚、1式の代わりにを用いることもできるが、舵感度時定数Ts3nTsnであるため、以下、簡略化のために舵感度時定数Ts3nを0とする。上式を状態空間表現で表すとになる。ここでT:転置,rは旋回角速度図1A参照をΔa、Δbはそれぞれパラメータ不確かさで、実際値とノミナル値との差異によってに定める。1.2.2外乱モデル外乱モデルは、図2に示すように、舵角オフセットモデルと波浪モデルとからなる。舵角オフセットモデルは風などに誘起された方位軸回りに作用するモーメントを舵角換算したものとする。波浪モデルは白色ノイズが入力した狭帯域フィルタの出力を方位換算したものとする。波浪モデルを状態空間表現で表すとになる。ここでδoは一定値とした舵角オフセット成分を、xw=ξ,ψwT,ξ:変数,ψwは平均値ゼロの有色の波浪成分を、νは白色ノイズN0,1をKwn,ζwn,ωwnはそれぞれ波浪モデルのゲイン、減衰係数と中心周波数とを示す。波浪モデルはパラメータ不確かさを含まないとする。1.3閉ループ制御系1.3.1フィードバック制御器フィードバック制御器12は、図2に示すように状態フィードバック制御器22と推定器オブザーバ24とからなる。閉ループ制御系は制御対象18とフィードバック制御器12とから構成しになる。ここで:推定値を示し、はフィードバックゲインをf1:比例ゲイン、f2:微分ゲイン、1:舵角オフセットに対応,は推定ゲインをを示し、Oはゼロ行列を示す。推定器出力はになる。1.3.2閉ループ特性行列閉ループ特性行列は、推定誤差をを用いて上式よりになる。ここでAcは特性行列を示す。船体モデルに状態フィードバックを行ったループをフィードバックループAc1,1、2次系と呼び,制御対象に推定ゲインのフィードバックを行ったループを推定ループAc2,2、5次系と呼ぶ。上式より要素Ac1,1のフィードバックループによる固有値と、要素Ac2,2の推定ループによる固有値は、ΔAs,ΔBsの影響を受ける。1.4ノミナル値による特性多項式フィードバックループ、推定ループおよび波浪モデルのノミナル値の特性多項式をとする。また、15式は、とする。ここでdet:行列式、s:ラプラス演算子、I:適当な単位行列、Df:フィードバックループの特性多項式添え字f、Dest:推定ループの特性多項式、Dw:波浪モデルの特性多項式添え字w、De、Dew、Deo:それぞれ船体モデル添え字e、波浪モデル添え字ewと舵角オフセットモデル添え字eoとに対応する特性多項式、ζ、ω:2次系のそれぞれ減衰係数、固有周波数を示す。1.5状態フィードバックゲイン状態フィードバック制御器22において、ノミナル値によるフィードバックループの状態フィードバックゲインを定める。状態フィードバックは図2に示すように、PD制御を用いており、とおき、設計パラメータにf1とζfとを選ぶと、f2とωfとはそれぞれとして得られる。1.6推定ゲイン2次の推定ゲインは、図3Aに示される外乱無しの船体モデルのみで設計される2次推定器におけるk1=k11、k2=k21である。3次の推定ゲインは、図3Bに示される舵角オフセットモデルを含めて設計される3次推定器におけるk1=k13、k2=k23、k5=k53である。4次の推定ゲインは、図3Cに示される波浪モデルを含めて設計される4次推定器におけるk1=k14、k2=k24、k3=k34、k4=k44である。5次の推定ゲインは、図3Dに示される波浪モデルと共に舵角オフセットモデルを含めて設計される5次推定器におけるk1=k15、k2=k25、k3=k35、k4=k45、k5=k55である。5次推定器の特性多項式は、になる。一方、設計する特性多項式は、15式を展開してになる。よって、推定ゲインは上記2つの多項式の係数から次式より求まる。ここでを示す。以降の設計に利用する推定ゲインを以下に示す。また、舵感度時定数Ts3nを考慮した場合には、以下となる。Desのパラメータはフィードバックループのωfを基準に、と定める。ここでρ>1を推定係数と呼ぶ。ρを閉ループ安定性に基づいて以降で設計する。1.7設計パラメータと制御パラメータ特性多項式の設計パラメータと制御パラメータを、以下表1にまとめる。上記表より、Desの減衰係数ζeと推定係数ρを与えれば、制御系が構成されることになる。尚、上記表から微分ゲインf2>0、ζf=1/√2と定めるため、23式からωfTsn>1/√2になる。船体パラメータのノミナル値Ksn、Tsnは、船体の旋回の度に同定器によって同定され、更新される。波浪モデルの減衰係数ζwnと固有角周波数ωwnも周波数同定器によって同定される。2設計前章で検討したように、設計は、推定器の減衰係数ζeと推定係数ρに帰着する。まず、特性根で推定根を不安定にさせるパラメータ不確かさを特定し、仕様を定める。次いで、減衰係数ζeを推定ループ単体から推定誤差が低減する適正値を設定する。そして閉ループからρに対する仕様付近の推定根の特性を明らかにし、解の予測値を設定し、ρの求解算法を推定根特性及び予測値から算出する。2.1パラメータ不確かさ2.1.1特性多項式の導出12式の特性行列の行列式である特性多項式は、ノミナル値による特性多項式とパラメータ不確かさによる多項式との和に整理される。特性多項式を因数分解すると特性根が得られる。すなわちになる。推定器の次数ii=0、2、3、4、5に対応した次数ciの特性多項式は次式になる。ここで、を示すi:次数。2.1.2パラメータ不確かさの影響パラメータ不確かさが閉ループ制御系に及ぼす影響を検討する。特性多項式Dcsにおいて、パラメータ不確かさΔa,Δbを開ループゲインと見なせば、根軌跡手法が適用でき、閉ループ安定性を検討できる。開ループ伝達関数GHsを31式よりに定めると、閉ループ特性は図4に示すものと等価になる。根軌跡手法を用いて、推定器の次数と漸近線の本数上式で分母の極数から分子のゼロ点数を差し引いた数との関係を表2にまとめる。表2より2次から5次までの推定器の漸近線の本数と形状は同じであり、推定器に起因する閉ループ特性は、2次推定器を用いたDc4sを基本とすることができることが分かるから、推定器による閉ループ特性を把握するためには、Dc4sを調べればよい。パラメータ不確かさΔa,Δb毎に対するDc4sの根軌跡をそれぞれ図5a及びbに示す。同図で×は極を、○はゼロ点を示す。この導出の条件として、ωfTsn=1.58Tsn=62.5、Ksn=0.04、f1=1とする。図5から次のことが分かる。Δa>0の場合はフィードバックループが不安定傾向になり、推定ループは間接的に影響する。よって、ζfの調整によって直接改善することができる。Δb>0の場合は推定ループが不安定傾向になり、推定根が直接影響し、フィードバックループは間接的に影響する。従って、推定係数ρの設定はΔbを利用することがよいことが分かる。2.1.3設計仕様推定係数ρが満足すべき仕様は、実際の制御対象の不確定要素、非線形性などを考慮して定める。具体的には、パラメータ不確かさの仕様Δspecは、Δbと対応するものとして、ある値第1仕様値に定める。そのため、パラメータ不確かさΔbの無次元量を導入する。ここで、’は無次元量を意味する。そして、例えば、パラメータ不確かさの仕様Δspecは、パラメータ不確かさΔbがKsn/Tsnの2倍に相当するものとして、に定める。また、閉ループ安定性の仕様ζspecは、特性根の中で最小の減衰係数に対応させて、ある値第2仕様値に決める。この値としては、例えば、0.4とすることができる。即ち、である。ここで上付き文字はパラメータ不確かさを含んだ場合を意味する。仕様ζspec=0.4はステップ応答のオーバシュート25に相当する量で、負の実軸に対して偏角cos-10.4=66.4[deg]になる図5b参照。2.2減数係数ζe外乱モデルを含まない2次推定器の減衰係数ζeはパラメータ不確かさに起因した推定誤差を低減させる値に選ぶ。船体モデルに対する推定ループは、9式及び図3Aを参照すると、になる。43式から次式を得る。上式及びそれを微分した式と実際値Ks、Tsの船体モデル1式参照とを44式に代入すると、次式になる。ここで、Δe2、Δe1はパラメータ不確かさを表す。上式から伝達関数は、になる。ここで、Desは、17式と等価である。また、48式を45式に代入すると、rの伝達関数は次式となる。伝達関数Gψψ^|2s、Gψr^|2sは、Δe2、Δe1によって共に推定誤差を生じるが、閉ループ安定性により強く関与するGψr^|2sを用いて減衰係数ζeを定める。2.2.1Δe1に関して50式と28式とより、Δe1に関連する項を抽出すると、となる。Δe1を微小として、ζee1をΔe1に関する級数で近似すると、となる。これより、ζee1がΔe1による誤差を受けないようにするため、gζe=0を解くと、次式を得る。2.2.2Δe2に関して同様に、50式と28式とより、Δe2に関連する項を抽出すると、となる。Δe2を微小として、ζee2をΔe2に関する級数で近似すると、となる。これより、ζee2がΔe2による誤差を受けないようにするための条件は、となる。2.2.3減数係数ζeの決定上記53式、56式に近い値を減数係数ζeの値として決定することが好ましい。具体的には、ωeTsnが1より十分大きいものと仮定すると、53式から減数係数ζeは、と定めることができ、この値は56式からも、さほど離れない値であると期待できる。2.3推定係数ρ推定係数ρは、31式のパラメータ不確かさ項に依存する。まずは、船体パラメータに対応する33式に基づいて、ρ2及び2次推定根の特性を把握する。そして、34~36式のパラメータ不確かさ項の係数を、33式のそれと比較することで、外乱モデルに起因するρii=3,4,5及び3~5次推定根の特性を把握する。手順としては、4次特性根とその2次推定根とを2次推定係数ρ2及びΔspecに関して解析し、ρ2と2次推定根との予測値を数値解法を利用して求める。5~7次の特性多項式は、4次の特性多項式に基づいて解析し、5~7次の特性根における推定係数と推定根との予測値をρ2と2次推定根とから求める。2.3.12次推定器31式の4次特性多項式は、になる。ここで、Dsは、パラメータ不確かさを含んだ多項式で、を示す。37式、23式及び27式より、を示す。パラメータ不確かさは、39式、22式を用いて、になり、ΔbB2sはf1が相殺され、f1の影響を受けないことが分かる。次に、以降の3次推定器の5次特性多項式の解析の準備として、2次推定根の特性について検討することにする。2次推定根は、図5bよりΔbの増加に伴って不安定側に移動し、その減衰係数ζe2は4次特性根の中で最小になる。その特性を把握するために、ωfTsn=1.58Tsn=62.5、Ksn=0.04、f1=1において、Δb’を1,2仕様値,3のそれぞれの場合に対して、ρ2を1~10まで増加させたときの数値解を求めて、ζf2、ωf2、ζe2、ωe2をそれぞれ求めると、図6に示すものとなる。この図6から、ρ2の増加に対応して、以下の傾向が把握される。ζe2は増加し、ζe=1/√2に漸近し、ωe2は単調増加する。ζf2、ωf2は徐々に変化しなくなり、一定値に近づく。Δb’が大きくなるに従って、ζe2、ωe2の線形範囲が拡大する。よって、を満足するρspec=ρ2は存在し図6の例では、ρ2=4.58、そのときωspec=ωe2とおく。図7は、仕様付近の2次推定根の特性を示す。Δb’/Δspec=1.0,1.1,1.2,1.3のそれぞれの場合の、ρ2’=ρ2/ρspecを1~1.4まで変化させたときの数値解を用いる。図7から以下のことが把握される。Δb’/Δspecが大きくなるほど、線形性は向上する。ζe2,ωe2はそれぞれのdζe2dρ,dωe2dρがほぼ一定と扱え、ζspec付近でほぼ直線近似が成り立つ。図8は図7の結果から導いたρ2およびζe2,ωe2の正規化した特性を示す。横軸にΔb’/Δspecを、縦軸にζe2’=ζe2/ζspec、ρ2’=ρ2/ρspec、ωe2’=ωe2/ωspecの無次元値を用いる。図8から以下のことが把握される。ζe2’はΔb’/Δspecにほぼ反比例する。このことから、ζe2’の予測値をと近似することができる。ここで、:予測値を表す。近似誤差=予測値÷計算値-1はΔb’/Δspec=1.3で約2になる。ρ2’、ωe2’はほぼ同値でΔb’/Δspecにほぼ比例する。図8の読取から、仕様値に対するパラメータ不確かさであるΔb’/Δspec、に対応するρ2’、ωe2’の予測値を次式で近似することができる。ここで近似誤差はΔb’/Δspec=1.3のときρ2’が約0.8、ωe2’が約1.8になる。このように、Δbの増加によるζe2の劣化、ωe2の変化およびρ2の増加は、Δb’/Δspecによる予測値から把握できることが分かる。2.3.23次推定器3次推定器は2次推定器に、舵角オフセット成分δoの推定を加えたもので、5次特性多項式は、34式から、になる。ここで、B3sは、28式より、を示す。2次推定器との対比を行って、2次推定器からの変化の特性を把握することによって、3次推定器の解析を行う。まず、65式のパラメータ不確かさの多項式を、の因数分解した形に置き換える。ここで、ζf=ζeを用いてを示す。2次推定器で用いた具体的数値とρ2=4.58、ρo=0.1から4azcz/bz2=0.305となり、4azcz/bz2は小さいことが推定される。よって、なる近似を用いると√1ε≒1ε/20<|ε|<1、Zsの根として次式を得るz2<z1とする。68式において分母は1に近いため、z1≒ωeoであり、しかも、Deo3sの特性根-ωeo3は極-ωeoからゼロ点z1の間にあるので、Deo3s≒Deosと近似することができる。これより65式は次式になる。71式を見ると、Dc43Δbsは2次推定器の58式のDc4Δbsと同じ形になっていることが分かる。そこで、両者のΔbの係数を比較するために、59式のs1係数b12、s0係数b02と71式のs1係数az、s0係数azz2とをそれぞれ比較すると、になる。72式及び73式の比がそれぞれ1より大きいことから、パラメータ不確かさΔbの係数は、3次推定器の5次特性多項式の方が、2次推定器の4次特性多項式に比べ、大きくなることが分かり、2次推定器におけるパラメータ不確かさΔbが等価的に大きくなったものとみなすことができる。そのため3次推定係数ρ3は、ρ2に対してパラメータ不確かさΔbの等価的増加分を修正する必要がある。Δbの係数比を具体的な数値を代入して計算するとρo=0.1になる。上式から大きい値のs1係数を選ぶこととすると、2次推定器の4次特性多項式においてパラメータ不確かさがazb12倍増加したことと等価になる。よって、前節の2次推定器の特性から64式を用いて、ρ’3の予測値を換算するとと定めることができる。これよりρ3、ωe3の予測値を次式に定めることができる。2.3.34次推定器4次推定器は2次推定器に波浪成分ψwの推定を加えたもので、6次特性多項式は35式から、になる。ここで、29式よりである。3次推定器での手法と同様に、2次推定器との対比を行って、4次推定器の解析を行う。そのため、B4の係数b14と、B2の係数b12との比較をすると、79式からとなる。詳細は省略するが、f1Tsnωe22/b12>1であるため、上記式は、kwに関して右下がりの直線となる。波浪パラメータの中心角周波数比をと定め、ζwn=0.1として数値計算するとを得る。最悪ケースのb14b12=1.121は74式のazb12=1.137とほぼ同等であるが、rw≫1でb14はb12に収束する。78式でDwとDewとは次の条件のときに近似できるので、2次推定器の特性と等価になり、4次推定係数ρ4≒ρ2,ωe4≒ωe2になる。上記の近似ができない場合、ρ4は減衰係数差ζe-ζwnに比例し、周波数比rwに反比例する傾向を持つとし、その予測値を次式に定める。2.3.45次推定器5次推定器は2次推定器に舵角オフセット成分δo、波浪成分ψwの推定を加えたもので、3次及び4次の推定根特性を踏襲するものとし、7次特性多項式は36式から、次式となる。上式の予測値を3次、4次の推定根特性から次式に定める。2.4波浪パラメータ条件Dc7Δbsにおいて、仕様を満たすがωe5=ρ5ωf>ωwnになった場合、波浪パラメータ条件表1に反する。そのときζwnζeに変更して再度ρ5を求める。ノッチフィルター効果外乱除去比ζwnζewを失うが、ρ5が小さくなるのでωe5を下げローパスフィルタ効果を得ることができる。すなわち85式でになり、Dc5Δbsの3次推定器の特性に変更される。2.5まとめ2.4節の結果をまとめると、推定根のζe2ωe2はζspec付近で、ρ2にほぼ比例する。外乱モデルを含むρ5とωe5との予測値は、外乱モデルを含まないρ2とωe2によって求められる。7次特性根からρ5を計算するとき、ρ5とζe5ωe5との初期値が適切であると求解が確実になる。以下、その求解算法について説明する。3.求解算法推定係数ρを求める算法を示す。求められたρはパラメータ不確かさを印加した閉ループ系において、仕様Δspecζspecを満足する。算法は確実な求解法、簡単な処理と少ない計算量とのため、前章で導き出した性質を用いて、3つの段階により構築する。3.1メインフロー次の流れを採用する図9158式でζe2=ζspecになるρ2,ωe2を求めるS10。285式でζe5=ζspecになるρ5,ωe5を求めるために、86式の予測値を初期値に設定するS20。またωe5>ωwnのときS30でyes、ζwn=ζeにしてS40、2を繰り返す。3求めたρ5をρ=ωe/ωfとして制御系を構成するS50。3.2ステップ2サブフロー1ρρ2、ρ5はその性質を利用して2段階で求める図10参照。1ζspecを挟む範囲のρを求める:ζeはρに対して単調増加になる性質を利用する。初期値はとおく。ステップ3に初期値を送るとζe,ωeが戻される。反復ループ添え字i:回数でを更新する。上記にてeζi=ζei-ζspec,eζi>0のときブレークし、下限値と上限値とを設定する図10参照。2次に、ζe=ζspecとなるρを求める。ζeはζspec付近でρにほぼ直線近似になる性質を利用する。反復ループ添え字i:回数でを更新する。上記にて|eζi|<εε:微小項のときブレークする。|eζi|≧εのときに、範囲を次式に置き換えて反複ループを実行する図10参照。3.3ステップ3初期値を用いて次の処理を行い、ζeとωeとを取り出す。14次式は代数解法によって、最小の減衰係数に対応する2次式の係数を求める。27次式は任意の数値解法によって求める。具体的には、ベアストウ法と組立除法とによって、2次式の係数を求め、残りの5次式を因数分解する場合はニュートン-ラプソン法初期値ωeoと代数解法とを用いる。なお、上記の解法は公知のため説明を省く。4.計算結果本発明を数値計算によって検証する。船体パラメータおよび制御パラメータを表3にまとめる。ただし設計パラメータについては表1を参照されたい。また、仕様は40、42式に従っている。4.12次推定器表3から以下のことが分かる。ρ2の変化は微分ゲインf2のそれより小さいが、ωfTsnに関係する。ρ2は仕様から判断すると実用範囲の値である。4.23次推定器図11を参照するとazb12:ωfTsn=1の場合は他の場合より15~30大きくなるが、ωfTsn>1の場合はほぼ近い値を示す。-azz2b02:近似値はazb12の90弱になり、ρoが大きくなると近似誤差が増加する。誤差ρ3’,誤差ωe3’:ρo≦0.2とすれば、両者とも5以内の誤差に収まる。よって、3次推定器を2次推定器に近似し予測値を求める方法は妥当であり、予測値ρ3’,ωe3’は有効であると言える。4.34次推定器図12ωfTsn=1.58を利用を参照するとb14b12、kw:共にrw=1で1に近く、rw=10でほぼ1に収斂する。ζe4’,ωe4’:B4sのみの場合波浪成分Dews,Dwsを省くは波浪成分を含む場合より変化量が小さい。図13より、ρ4’の予測値はζwn=0.1かつrwが低減で誤差を生じるが、それ以外では誤差は小さいと見なせる。ωe4’の予測値は一定のため、低い領域で誤差を生じる。これに対して、波浪成分による影響が支配的になり、予測値ρ4’,ωe4’は実用精度をもつ。4.45次推定器図14は図13にρoを追加した場合で、ほぼ同様の傾向を示し、オフセットは予測値とほぼ等価である。図15はζwn=0.1,ωfTsnを変化させた場合で、ωfTsn=1のrwの低減で急増する。予測値ρ5’,ωe5’はrwの低減で波浪成分誤差をもつが、波浪パラメータ条件破線より右側部分が条件を満足する領域から予測値として実用精度をもつことが分かる。4.5計算時間の比較本発明の求解算法の計算時間を標準算法と比較する。標準算法は既知の推定係数ρsolから探索範囲[ρlow,ρhigh]=ρsol[1cp,cp]を与えてMatlabのfminbndを利用、12式の特性行列から特性根を求め、最小の減衰係数を取り出し、それが仕様と一致する推定係数を求める。計算結果を表4にまとめる。条件:図14,ζwn=0.1である。本発明の計算時間は0.582,0.581,0.585s,平均0.583sになり、標準算法の計算時間に比べ、1桁程度速いことが確認できた。5.まとめ推定器を含んだ保針ループに対して、推定根に影響を与えるパラメータ不確かさを特定し、仕様を定め、推定根特性を仕様付近で把握することで、ρの予測値を求めた。外乱モデルなしの2次推定根特性を基に外乱モデルありの推定根特性を解析した。求解算法は数値解と推定根特性とを利用し予測値を初期値に採用した。数値計算によって、推定根特性と予測値との妥当性を確認し、求解算法の有効性を標準算法の計算時間と比べ確認した。"