船舶用自動操舵装置
【課題】パラメータの不確かさの存在下で閉ループの安定性を図ることができる船舶用自動操舵装置を提供する。
【解決手段】制御対象１８からの検出方位から外乱を除去した推定値を出力する推定器２４と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器２２とを有するフィードバック制御器１２を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象１８と共に閉ループを構成する。フィードバックループの特性多項式と、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなる推定器２４の特性多項式からなる閉ループの特性多項式により、推定器の船体モデルの固有角周波数ωeを設定する。
【解決手段】制御対象１８からの検出方位から外乱を除去した推定値を出力する推定器２４と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器２２とを有するフィードバック制御器１２を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象１８と共に閉ループを構成する。フィードバックループの特性多項式と、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなる推定器２４の特性多項式からなる閉ループの特性多項式により、推定器の船体モデルの固有角周波数ωeを設定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、推定器を有する船舶用自動操舵装置に関し、特にパラメータの不確かさの存在下で閉ループの安定性を図ることができる船舶用自動操舵装置に関する。
【背景技術】
【０００２】
船舶用自動操舵装置は、一定の設定方位に、ジャイロコンパスから検出される船首方位を追従させるために舵を制御する装置であり、保針と変針との機能を持つ。
【０００３】
その保針制御系は、大きな慣性力の船体を小さな制御力の操舵機によって制御し、かつ検出方位に含まれる外乱の大きさが方位と同等かそれ以上になる特徴をもつため、受動的な制御を実現するフィードバック制御器で制御される。図２に示したように、フィードバック制御器１２は、小さな比例ゲインを基にしたＰＤ制御を行う状態フィードバック制御器２２と、状態推定と外乱除去とを行う推定器２４から構成される。閉ループ安定性と外乱除去性への影響及び設計自由度は、実用上状態フィードバック制御器２２よりも推定器２４の方が大きい。
【０００４】
制御対象１８の船体モデルは、載荷の変化によるノミナル値のパラメータ不確かさを持ち、そのためノミナル値のモデルベースで構成する推定器２４は載荷の変化によって推定値に誤差を生じる。このため、推定値の状態フィードバックによる閉ループ安定性は、推定誤差によって劣化し、船体を蛇行航行させる（ヨーイング）おそれがある。
【０００５】
本発明者らは、かかる問題を解決するために、特許文献１において、船体パラメータのノミナル値のパラメータ不確かさを積極的に考慮に入れて、該パラメータ不確かさを起因とする推定誤差を小さくすることができるようにした船舶用自動操舵装置及びその推定器の設計方法を提案している。
ここでは、推定器の特性多項式として、
【０００６】
【数１】
と置き、λe2は、波浪モデルが仮に無いとしたときの船体モデルの状態量を推定するための特性多項式であり、ζe、ωeがそれぞれ船体モデルの状態量を推定するための減衰係数、固有周波数としたときに、ωeを操舵系周波数ωfのρ（＞１）倍(以下、ρ：推定係数)に設定し、ζeを１／√２に設定するようにして推定器を設計することにより、閉ループ安定性を低下させる要因となる推定角速度のパラメータ不確かさによる推定誤差を低減させることができる、ことを見出した。
【０００７】
さらに、特許文献１では、船体モデルと２つの外乱モデル（波浪モデル、舵角オフセットモデル）とからなる制御対象に対して、推定誤差が船体モデルのみのときの推定誤差と等価になるように推定器の固有周波数を外乱モデル特性から修正することを、提案している。これによれば、推定器の制御対象は船体モデルのみとなり、外乱モデルを無視することができるようになる。
【０００８】
また、特許文献２では、推定係数ρ＝ωe／ωfを、操舵ループの代表根に着目し、代表根の減衰係数とパラメータ不確かさとの関係から制御ゲインを設計している。即ち、ここでは推定器がない閉ループ系において、フィードバックゲインの設定、操舵系の特性およびパラメータ不確かさの影響を調べ、推定器がある閉ループ系において、特性多項式の導出、代表根の移動、代表根の減衰係数を調べて数値計算を実施した結果、パラメータ不確かさの影響を見積もることによって、推定器と操舵系との固有周波数を関係づける推定係数ρを設定している。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特許第４６０３８３２号公報
【特許文献２】特開２００７−１１８８２８号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかしながら、従来においては、パラメータ不確かさおよび外乱モデルの影響を間接的、近似的に扱うために、閉ループ安定性に多少の変動を持つという問題がある。即ち、推定ループにおいて閉ループ安定性を劣化させる場合があるが、従来の設計による装置では、推定ループにおけるパラメータ不確かさの影響を直接扱っていない。
【００１１】
また、船体モデルはパラメータ不確かさを複数もつために、特許文献２で採用したパラメータ不確かさ以外が船種、海象、載貨などにより無視できない状況を生じる場合があり、用いたパラメータ不確かさによる影響が減衰係数に間接的に作用するために局所的な安定性しか確保することができない場合がある、という問題がある。
【００１２】
本発明はかかる課題に鑑みなされたもので、外乱モデルを含んだ推定器と操舵系の閉ループにおいて、パラメータ不確かさに起因する推定誤差を低減することができる推定ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。
【００１３】
また、本発明の更なる目的は、推定ループの根を不安定にさせるパラメータ不確かさを特定し、それを用いて制御ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。
【課題を解決するための手段】
【００１４】
前述した目的を達成するために、請求項１記載の発明は、設定方位ψＲに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、
【００１５】
【数２】
前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、
【００１６】
【数３】
船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔを、以下で表したときに、
【００１７】
【数４】
【００１８】
パラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする。
【００１９】
請求項２記載の発明は、請求項１記載の船舶用自動操舵装置において、前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄc4Δを、
【００２０】
【数５】
としたときに、Ｄc4Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζe2＊とし、ζe2＊が第２仕様値（ζspec）の値を満足するωe2／ωf（＝ρ2）を求め、このρ2を用いた初期値とした数値解法によって前記ωe／ωf（＝ρ）を求めることを特徴とする。
【００２１】
請求項３記載の発明は、請求項２記載の船舶用自動操舵装置において、前記初期値は、
【００２２】
【数６】
で表されることを特徴とする。
【００２３】
請求項４記載の発明は、請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置において、前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさであることを特徴とする。
【００２４】
請求項５記載の発明は、請求項４記載の船舶用自動操舵装置において、前記特性多項式ＤcΔは、
【００２５】
【数７】
で表され、ここで、
Ｄw（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、
【００２６】
【数８】
と表され、Ｂ（ｓ）は、
【００２７】
【数９】
で表されることを特徴とする。
【発明の効果】
【００２８】
本発明によれば、外乱モデルとしての波浪モデル及び舵角オフセットモデルを含んだ推定器を含む閉ループの特性多項式に基づき、特性多項式に船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔからパラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することから、閉ループ安定性を増加させることができる。
【００２９】
また、前記パラメータ不確かさΔを、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさとすることで、閉ループ安定性に及ぼす影響の高いパラメータ不確かさに起因する推定根の不安定性を排除し、従来のものより信頼性を向上させることができる。
【００３０】
外乱モデルを含まない推定器の推定根から出発して、外乱モデルを含んだ推定器の推定根の予測値を求めて、特性多項式ＤcΔの特性根を求めることで、求解を高い収束性で短時間に行うことができるようになる。
【図面の簡単な説明】
【００３１】
【図１】本発明による船舶用自動操舵装置及び制御対象の全体構成を表すブロック図である。
【図１Ａ】船体モデルのブロック図である。ただし、舵感度時定数Ｔs3n項を省く。
【図２】閉ループ制御系を表すブロック図である。
【図３Ａ】２次推定器のブロック図である。
【図３Ｂ】３次推定器のブロック図である。
【図３Ｃ】４次推定器のブロック図である。
【図３Ｄ】５次推定器のブロック図である。
【図４】開ループ伝達関数による閉ループを表すブロック図である。
【図５】Ｄc4(s)の根軌跡を表すｓ平面である。
【図６】４次特性多項式の数値解から求めたρ2に対するζf2*、ωf2*、ζe2*、ωe2*とこれらの変化率の特性を表すグラフである。
【図７】仕様値付近での正規化したρ2に対するζe2*、ωe2*とこれらの変化率の特性を表すグラフである。
【図８】図６の結果から導いたρ2、ζe2*、ωe2*の正規化した場合の特性を表すグラフである。
【図９】本発明による求解法を表すフローチャートである。
【図１０】ρを求める原理を表す図である。
【図１１】３次推定器の推定根の特性を表す数値結果である。
【図１２】４次推定器の推定根の特性を表す数値結果である。
【図１３】４次推定器の推定根の特性を表す数値結果である。
【図１４】５次推定器の推定根の特性を表す数値結果である。
【図１５】５次推定器の推定根の特性を表す数値結果である。
【発明を実施するための形態】
【００３２】
以下、図面を用いて本発明の実施の形態を説明する。
【００３３】
１．１ 保針系の構成
本発明で対象とする船舶用自動操舵装置１０の保針モードは、図１に示すように操舵機１４及び船体１６を含む制御対象１８を設定方位ψRに船首方位ψを追従させるため、フィードバック制御器１２で偏差
【００３４】
【数１０】
から演算された舵角指令δcによって制御する。δcは操舵機１４の舵角δを経て船体に作用する。船体はパラメータ不確かさをもち外乱が印加する。操舵機１４は舵角および舵速度の振幅制限をもつが、保針時のδcおよび
の振幅が制限以下と見なせるので操舵機を省略する（δ＝δc）。なお設定方位ψRは保針時一定になるため、簡単化のためψR＝０とする。
【００３５】
１．２ 制御対象とパラメータ不確かさ
図１に示した制御対象１８は、操舵機１４を除き、船体モデルと外乱モデルとから構成される。
【００３６】
１．２．１ 船体モデル
船体モデルは、図１Ａに示される野本の線形モデル
【００３７】
【数１１】
を用いる。ここで、ｓはラプラス演算子を、ψは船首方位を、ｒを旋回角速度を、Ｋsn，Ｔsn，は操縦性指数（船体パラメータ）でそれぞれ旋回力ゲイン、時定数のノミナル値（添え字ｓ：船体、ｎ：ノミナル値）を示す。尚、（１）式の代わりに
【００３８】
【数１２】
を用いることもできるが、舵感度時定数Ｔs3n<Ｔsnであるため、以下、簡略化のために舵感度時定数Ｔs3nを０とする。
【００３９】
上式を状態空間表現で表すと
【００４０】
【数１３】
になる。ここで
【００４１】
（Ｔ：転置），ｒは旋回角速度（図１Ａ参照）を
【００４２】
【数１４】
【００４３】
【数１５】
Δa、Δbはそれぞれパラメータ不確かさで、実際値とノミナル値との差異によって
【００４４】
【数１６】
に定める。
【００４５】
１．２．２ 外乱モデル
外乱モデルは、図２に示すように、舵角オフセットモデルと波浪モデルとからなる。舵角オフセットモデルは風などに誘起された方位軸回りに作用するモーメントを舵角換算したものとする。波浪モデルは白色ノイズが入力した狭帯域フィルタの出力を方位換算したものとする。
【００４６】
波浪モデルを状態空間表現で表すと
【００４７】
【数１７】
【００４８】
【数１８】
になる。ここでδoは一定値とした舵角オフセット成分を、ｘw ＝ [ξ，ψw ]T ，ξ：変数，ψw は平均値ゼロの有色の波浪成分を、νは白色ノイズＮ（０，１）を
【００４９】
【数１９】
Ｋwn，ζwn，ωwn はそれぞれ波浪モデルのゲイン、減衰係数と中心周波数とを示す。波浪モデルはパラメータ不確かさを含まないとする。
【００５０】
１．３ 閉ループ制御系
１．３．１ フィードバック制御器
フィードバック制御器１２は、図２に示すように状態フィードバック制御器２２と推定器（オブザーバ）２４とからなる。閉ループ制御系は制御対象１８とフィードバック制御器１２とから構成し
【００５１】
【数２０】
になる。ここで
【００５２】
【数２１】
（^：推定値）を示し、
はフィードバックゲインを（ｆ1：比例ゲイン、ｆ2：微分ゲイン、１：舵角オフセットに対応），
は推定ゲインを
【００５３】
【数２２】
を示し、Ｏはゼロ行列を示す。推定器出力は
【００５４】
【数２３】
になる。
【００５５】
１．３．２ 閉ループ特性行列
閉ループ特性行列は、推定誤差を
【００５６】
【数２４】
を用いて上式より
【００５７】
【数２５】
【００５８】
【数２６】
になる。ここでＡcは特性行列を示す。船体モデルに状態フィードバックを行ったループをフィードバックループ（Ａc（１，１）、２次系）と呼び，制御対象に推定ゲインのフィードバックを行ったループを推定ループ（Ａc（２，２）、５次系）と呼ぶ。上式より要素Ａc（１，１）のフィードバックループによる固有値と、要素Ａc（２，２）の推定ループによる固有値は、ΔＡs，ΔＢsの影響を受ける。
【００５９】
１．４ ノミナル値による特性多項式
フィードバックループ、推定ループおよび波浪モデルのノミナル値の特性多項式を
【００６０】
【数２７】
【００６１】
【数２８】
【００６２】
【数２９】
とする。また、（１５）式は、
【００６３】
【数３０】
【００６４】
【数３１】
【００６５】
【数３２】
とする。ここでｄｅｔ：行列式、ｓ：ラプラス演算子、Ｉ：適当な単位行列、Ｄf：フィードバックループの特性多項式（添え字ｆ）、Ｄest：推定ループの特性多項式、Ｄw：波浪モデルの特性多項式（添え字ｗ）、Ｄe、Ｄew、Ｄeo：それぞれ船体モデル（添え字ｅ）、波浪モデル（添え字ｅｗ）と舵角オフセットモデル（添え字ｅｏ）とに対応する特性多項式、ζ、ω：２次系のそれぞれ減衰係数、固有周波数を示す。
【００６６】
１．５ 状態フィードバックゲイン
状態フィードバック制御器２２において、ノミナル値によるフィードバックループの状態フィードバックゲインを定める。状態フィードバックは図２に示すように、ＰＤ制御を用いており、
【００６７】
【数３３】
とおき、設計パラメータにｆ1とζfとを選ぶと、ｆ2とωfとはそれぞれ
【００６８】
【数３４】
【００６９】
【数３５】
として得られる。
【００７０】
１．６ 推定ゲイン
２次の推定ゲインは、図３Ａに示される外乱無しの船体モデルのみで設計される２次推定器におけるｋ1（＝ｋ11）、ｋ2（＝ｋ21）である。
【００７１】
３次の推定ゲインは、図３Ｂに示される舵角オフセットモデルを含めて設計される３次推定器におけるｋ1（＝ｋ13）、ｋ2（＝ｋ23）、ｋ5（＝ｋ53）である。
【００７２】
４次の推定ゲインは、図３Ｃに示される波浪モデルを含めて設計される４次推定器におけるｋ1（＝ｋ14）、ｋ2（＝ｋ24）、ｋ3（＝ｋ34）、ｋ4（＝ｋ４4）である。
【００７３】
５次の推定ゲインは、図３Ｄに示される波浪モデルと共に舵角オフセットモデルを含めて設計される５次推定器におけるｋ1（＝ｋ15）、ｋ2（＝ｋ25）、ｋ3（＝ｋ35）、ｋ4（＝ｋ45）、ｋ5（＝ｋ55）である。
【００７４】
５次推定器の特性多項式は、
【００７５】
【数３６】
になる。一方、設計する特性多項式は、（１５）式を展開して
【００７６】
【数３７】
になる。
よって、推定ゲインは上記２つの多項式の係数から次式より求まる。
【００７７】
【数３８】
ここで
【００７８】
【数３９】
を示す。
【００７９】
以降の設計に利用する推定ゲインを以下に示す。
【００８０】
【数４０】
【００８１】
【数４１】
【００８２】
【数４２】
【００８３】
また、舵感度時定数Ｔs3nを考慮した場合には、以下となる。
【００８４】
【数４３】
【００８５】
【数４４】
【００８６】
【数４５】
【００８７】
Ｄe（ｓ）のパラメータはフィードバックループのωfを基準に、
【００８８】
【数４６】
と定める。ここでρ＞１を推定係数と呼ぶ。ρを閉ループ安定性に基づいて以降で設計する。
【００８９】
１．７ 設計パラメータと制御パラメータ
特性多項式の設計パラメータと制御パラメータを、以下表１にまとめる。
【００９０】
【表１】
【００９１】
上記表より、Ｄe（ｓ）の減衰係数ζeと推定係数ρを与えれば、制御系が構成されることになる。尚、上記表から微分ゲインｆ2＞０、ζf＝１／√２と定めるため、（２３）式からωfＴsn＞１／√２になる。船体パラメータのノミナル値Ｋsn、Ｔsnは、船体の旋回の度に同定器によって同定され、更新される。波浪モデルの減衰係数ζwnと固有角周波数ωwnも周波数同定器によって同定される。
【００９２】
２ 設計
前章で検討したように、設計は、推定器の減衰係数ζeと推定係数ρに帰着する。まず、特性根で推定根を不安定にさせるパラメータ不確かさを特定し、仕様を定める。次いで、減衰係数ζeを推定ループ単体から推定誤差が低減する適正値を設定する。そして閉ループからρに対する仕様付近の推定根の特性を明らかにし、解の予測値を設定し、ρの求解算法を推定根特性及び予測値から算出する。
【００９３】
２．１ パラメータ不確かさ
２．１．１ 特性多項式の導出
（１２）式の特性行列の行列式である特性多項式は、ノミナル値による特性多項式とパラメータ不確かさによる多項式との和に整理される。特性多項式を因数分解すると特性根が得られる。すなわち
【００９４】
【数４７】
になる。推定器の次数ｉ（ｉ＝０、２、３、４、５）に対応した次数ｃｉの特性多項式は次式になる。
【００９５】
【数４８】
【００９６】
【数４９】
【００９７】
【数５０】
【００９８】
【数５１】
【００９９】
【数５２】
ここで、
【０１００】
【数５３】
を示す（ｉ：次数）。
【０１０１】
２．１．２ パラメータ不確かさの影響
パラメータ不確かさが閉ループ制御系に及ぼす影響を検討する。特性多項式Ｄc（ｓ）において、パラメータ不確かさΔa，Δbを開ループゲインと見なせば、根軌跡手法が適用でき、閉ループ安定性を検討できる。開ループ伝達関数ＧＨ（ｓ）を（３１）式より
【０１０２】
【数５４】
に定めると、閉ループ特性は図４に示すものと等価になる。根軌跡手法を用いて、推定器の次数と漸近線の本数（上式で分母の極数から分子のゼロ点数を差し引いた数）との関係を表２にまとめる。
【０１０３】
【表２】
【０１０４】
表２より２次から５次までの推定器の漸近線の本数と形状は同じであり、推定器に起因する閉ループ特性は、２次推定器を用いたＤc4（ｓ）を基本とすることができることが分かるから、推定器による閉ループ特性を把握するためには、Ｄc4（ｓ）を調べればよい。
【０１０５】
パラメータ不確かさΔa，Δb毎に対するＤc4（ｓ）の根軌跡をそれぞれ図５（ａ）及び（ｂ）に示す。同図で×は極を、○はゼロ点を示す。この導出の条件として、ωfＴsn ＝１．５８（Ｔsn ＝６２．５、Ｋsn ＝０．０４、ｆ1＝１）とする。
【０１０６】
図５から次のことが分かる。
・Δa＞０の場合はフィードバックループが不安定傾向になり、推定ループは間接的に影響する。よって、ζfの調整によって直接改善することができる。
・Δb＞０の場合は推定ループが不安定傾向になり、推定根が直接影響し、フィードバックループは間接的に影響する。
【０１０７】
従って、推定係数ρの設定はΔbを利用することがよいことが分かる。
【０１０８】
２．１．３ 設計仕様
推定係数ρが満足すべき仕様は、実際の制御対象の不確定要素、非線形性などを考慮して定める。
【０１０９】
具体的には、パラメータ不確かさの仕様Δspecは、Δbと対応するものとして、ある値（第１仕様値）に定める。そのため、パラメータ不確かさΔbの無次元量を導入する。
【０１１０】
【数５５】
ここで、（・）’は無次元量を意味する。
【０１１１】
そして、例えば、パラメータ不確かさの仕様Δspecは、パラメータ不確かさΔbがＫsn／Ｔsnの２倍に相当するものとして、
【０１１２】
【数５６】
に定める。
【０１１３】
また、閉ループ安定性の仕様ζspecは、特性根の中で最小の減衰係数に対応させて、ある値（第２仕様値）に決める。この値としては、例えば、０．４とすることができる。即ち、
【０１１４】
【数５７】
【０１１５】
【数５８】
である。ここで上付き文字(*) はパラメータ不確かさを含んだ場合を意味する。仕様ζspec＝０．４はステップ応答のオーバシュート２５%に相当する量で、負の実軸に対して偏角cos−１（０．４）＝６６．４［ｄｅｇ］になる（図５（ｂ）参照）。
【０１１６】
２．２ 減数係数ζe
外乱モデルを含まない２次推定器の減衰係数ζeはパラメータ不確かさに起因した推定誤差を低減させる値に選ぶ。
【０１１７】
船体モデルに対する推定ループは、（９）式及び図３Ａを参照すると、
【０１１８】
【数５９】
【０１１９】
【数６０】
になる。（４３）式から次式を得る。
【０１２０】
【数６１】
【０１２１】
上式及びそれを微分した式と実際値Ｋs、Ｔsの船体モデル（（１）式参照）
【０１２２】
【数６２】
とを（４４）式に代入すると、次式になる。
【０１２３】
【数６３】
ここで、Δe2、Δe1はパラメータ不確かさを表す。
【０１２４】
上式から伝達関数は、
【０１２５】
【数６４】
【０１２６】
【数６５】
になる。ここで、Ｄe（ｓ）は、（１７）式と等価である。
【０１２７】
また、（４８）式を（４５）式に代入すると、ｒの伝達関数は次式となる。
【０１２８】
【数６６】
【０１２９】
伝達関数Ｇψψ＾｜2（ｓ）、Ｇψｒ＾｜2（ｓ）は、Δe2、Δe1によって共に推定誤差を生じるが、閉ループ安定性により強く関与するＧψｒ＾｜2（ｓ）を用いて減衰係数ζeを定める。
【０１３０】
２．２．１ Δe1に関して
（５０）式と（２８）式とより、Δe1に関連する項を抽出すると、
【０１３１】
【数６７】
となる。Δe1を微小として、ζe*e1をΔe1に関する級数で近似すると、
【０１３２】
【数６８】
となる。これより、ζe*e1がΔe1による誤差を受けないようにするため、ｇ（ζe）＝０を解くと、次式を得る。
【０１３３】
【数６９】
【０１３４】
２．２．２ Δe2に関して
同様に、（５０）式と（２８）式とより、Δe2に関連する項を抽出すると、
【０１３５】
【数７０】
となる。Δe2を微小として、ζe*e2をΔe2に関する級数で近似すると、
【０１３６】
【数７１】
となる。これより、ζe*e2がΔe2による誤差を受けないようにするための条件は、
【０１３７】
【数７２】
となる。
【０１３８】
２．２．３ 減数係数ζeの決定
上記（５３）式、（５６）式に近い値を減数係数ζeの値として決定することが好ましい。具体的には、ωeＴsnが１より十分大きいものと仮定すると、（５３）式から減数係数ζeは、
【０１３９】
【数７３】
と定めることができ、この値は（５６）式からも、さほど離れない値であると期待できる。
【０１４０】
２．３ 推定係数ρ
推定係数ρは、（３１）式のパラメータ不確かさ項に依存する。まずは、船体パラメータに対応する（３３）式に基づいて、ρ2及び２次推定根の特性を把握する。そして、（３４）〜（３６）式のパラメータ不確かさ項の係数を、（３３）式のそれと比較することで、外乱モデルに起因するρｉ（ｉ＝３，４，５）及び３〜５次推定根の特性を把握する。
【０１４１】
手順としては、
・４次特性根とその２次推定根とを２次推定係数ρ2及びΔspecに関して解析し、ρ2と２次推定根との予測値を数値解法を利用して求める。
・５〜７次の特性多項式は、４次の特性多項式に基づいて解析し、５〜７次の特性根における推定係数と推定根との予測値をρ2と２次推定根とから求める。
【０１４２】
２．３．１ ２次推定器
（３１）式の４次特性多項式は、
【０１４３】
【数７４】
になる。ここで、Ｄ*（ｓ）は、パラメータ不確かさを含んだ多項式で、
【０１４４】
【数７５】
を示す。（３７）式、（２３）式及び（２７）式より、
【０１４５】
【数７６】
を示す。
【０１４６】
パラメータ不確かさは、（３９）式、（２２）式を用いて、
【０１４７】
【数７７】
になり、ΔbＢ2（ｓ）はｆ1が相殺され、ｆ1の影響を受けないことが分かる。
【０１４８】
次に、以降の３次推定器の５次特性多項式の解析の準備として、２次推定根の特性について検討することにする。２次推定根は、図５（ｂ）よりΔbの増加に伴って不安定側に移動し、その減衰係数ζe2*は４次特性根の中で最小になる。その特性を把握するために、ωfＴsn ＝１．５８（Ｔsn ＝６２．５、Ｋsn ＝０．０４、ｆ1＝１）において、Δb’を１，２（仕様値），３のそれぞれの場合に対して、ρ2を１〜１０まで増加させたときの数値解を求めて、ζf2*、ωf2*、ζe2*、ωe2*をそれぞれ求めると、図６に示すものとなる。この図６から、ρ2の増加に対応して、以下の傾向が把握される。
・ζe2*は増加し、ζe＝１／√２に漸近し、ωe2*は単調増加する。
・ζf2*、ωf2*は徐々に変化しなくなり、一定値に近づく。
・Δb’が大きくなるに従って、ζe2*、ωe2*の線形範囲が拡大する。
【０１４９】
よって、
を満足するρspec＝ρ2は存在し（図６の例では、ρ2＝４．５８）、そのときωspec＝ωe2*とおく。
【０１５０】
図７は、仕様付近の２次推定根の特性を示す。Δb’／Δspec＝１．０，１．１，１．２，１．３のそれぞれの場合の、ρ2’＝ρ2／ρspecを１〜１．４まで変化させたときの数値解を用いる。図７から以下のことが把握される。
・Δb’／Δspecが大きくなるほど、線形性は向上する。
・ζe2*，ωe2*はそれぞれのｄζe2*/ｄρ，ｄωe2*/ｄρがほぼ一定と扱え、ζspec付近でほぼ直線近似が成り立つ。
【０１５１】
【数７８】
【０１５２】
【数７９】
【０１５３】
図８は図７の結果から導いたρ2およびζe2*，ωe2*の正規化した特性を示す。横軸にΔb’ ／Δspecを、縦軸にζe2*’＝ζe2*／ζspec、ρ2’＝ρ2／ρspec、ωe2*’＝ωe2*／ωspecの無次元値を用いる。図８から以下のことが把握される。
・ζe2*’はΔb’ ／Δspecにほぼ反比例する。このことから、ζe2*’の予測値を
【０１５４】
【数８０】
と近似することができる。ここで、~：予測値を表す。
【０１５５】
近似誤差（＝（予測値）÷（計算値）−１）はΔb’／Δspec＝１．３で約２％になる。
・ρ2’、ωe2*’はほぼ同値でΔb’／Δspecにほぼ比例する。
【０１５６】
図８の読取から、仕様値に対するパラメータ不確かさであるΔb’ ／Δspec、に対応するρ2’、ωe2*’の予測値を次式で近似することができる。
【０１５７】
【数８１】
ここで近似誤差はΔb’／Δspec＝１．３のときρ~2’が約０．８％、ω~e2*’が約１．８％になる。
【０１５８】
このように、Δbの増加によるζe2*の劣化、ωe2*の変化およびρ2の増加は、Δb’／Δspecによる予測値から把握できることが分かる。
【０１５９】
２．３．２ ３次推定器
３次推定器は２次推定器に、舵角オフセット成分δoの推定を加えたもので、５次特性多項式は、（３４）式から、
【０１６０】
【数８２】
になる。ここで、Ｂ3（ｓ）は、（２８）式より、
【０１６１】
【数８３】
を示す。
【０１６２】
２次推定器との対比を行って、２次推定器からの変化の特性を把握することによって、３次推定器の解析を行う。まず、（６５）式のパラメータ不確かさの多項式を、
【０１６３】
【数８４】
の因数分解した形に置き換える。ここで、ζf ＝ ζe を用いて
【０１６４】
【数８５】
を示す。２次推定器で用いた具体的数値とρ2＝４．５８、ρo＝０．１から４ａzｃz／ｂz2＝０．３０５となり、４ａzｃz／ｂz2は小さいことが推定される。よって、
【０１６５】
【数８６】
なる近似を用いると（√（１＋ε）≒１＋ε／２ （０＜｜ε｜＜１）））、Ｚ（ｓ）の根として次式を得る（ｚ2＜ｚ1とする）。
【０１６６】
【数８７】
【０１６７】
【数８８】
【０１６８】
（６８）式において分母は１に近いため、ｚ1≒ωeoであり、しかも、Ｄeo3＊（ｓ）の特性根−ωeo3*は極−ωeoからゼロ点ｚ1の間にあるので、Ｄeo3*（ｓ）≒Ｄeo（ｓ）と近似することができる。これより（６５）式は次式になる。
【０１６９】
【数８９】
【０１７０】
【数９０】
【０１７１】
（７１）式を見ると、Ｄc43Δb（ｓ）は２次推定器の（５８）式のＤc4Δb（ｓ）と同じ形になっていることが分かる。そこで、両者のΔbの係数を比較するために、（５９）式のｓ1 係数ｂ12、ｓ0 係数ｂ02と（７１）式のｓ1 係数ａz、ｓ0 係数ａzｚ2とをそれぞれ比較すると、
【０１７２】
【数９１】
【０１７３】
【数９２】
【０１７４】
になる。（７２）式及び（７３）式の比がそれぞれ１より大きいことから、パラメータ不確かさΔbの係数は、３次推定器の５次特性多項式の方が、２次推定器の４次特性多項式に比べ、大きくなることが分かり、２次推定器におけるパラメータ不確かさΔbが等価的に大きくなったものとみなすことができる。そのため３次推定係数ρ3は、ρ2に対してパラメータ不確かさΔbの等価的増加分を修正する必要がある。Δbの係数比を具体的な数値を代入して計算すると（ρo＝０．１）
【０１７５】
【数９３】
になる。上式から大きい値のｓ1 係数を選ぶこととすると、２次推定器の４次特性多項式においてパラメータ不確かさがａz /ｂ12倍増加したことと等価になる。よって、前節の２次推定器の特性から（６４）式を用いて、ρ’3の予測値を換算すると
【０１７６】
【数９４】
と定めることができる。これよりρ3、ωe3*の予測値を次式に定めることができる。
【０１７７】
【数９５】
【０１７８】
２．３．３ ４次推定器
４次推定器は２次推定器に波浪成分ψwの推定を加えたもので、６次特性多項式は（３５）式から、
【０１７９】
【数９６】
になる。ここで、（２９）式より
【０１８０】
【数９７】
である。
【０１８１】
３次推定器での手法と同様に、２次推定器との対比を行って、４次推定器の解析を行う。そのため、Ｂ4の係数ｂ14と、Ｂ2の係数ｂ12との比較をすると、（７９）式から
【０１８２】
【数９８】
となる。詳細は省略するが、ｆ1Ｔsnωe22／ｂ12＞１であるため、上記式は、ｋwに関して右下がりの直線となる。
【０１８３】
波浪パラメータの中心角周波数比を
【０１８４】
【数９９】
と定め、ζwn＝０．１として数値計算すると
を得る。最悪ケースのｂ14 /ｂ12＝１．１２１は（７４）式のａz /ｂ12＝１．１３７とほぼ同等であるが、ｒw≫１でｂ14はｂ12に収束する。
【０１８５】
（７８）式でＤw(s)とＤew(s)とは次の条件のとき
【０１８６】
【数１００】
に近似できるので、２次推定器の特性と等価になり、４次推定係数ρ4≒ρ2，ωe4*≒ωe2*になる。上記の近似ができない場合、ρ4は減衰係数差（ζe−ζwn）に比例し、周波数比ｒwに反比例する傾向を持つとし、その予測値を次式に定める。
【０１８７】
【数１０１】
【０１８８】
２．３．４ ５次推定器
５次推定器は２次推定器に舵角オフセット成分δo、波浪成分ψwの推定を加えたもので、３次及び４次の推定根特性を踏襲するものとし、７次特性多項式は（３６）式から、次式となる。
【０１８９】
【数１０２】
【０１９０】
上式の予測値を３次、４次の推定根特性から次式に定める。
【０１９１】
【数１０３】
【０１９２】
２．４ 波浪パラメータ条件
Ｄc7Δb（ｓ）において、仕様を満たすがωe5 ＝ ρ5ωf＞ωwn になった場合、波浪パラメータ条件（表１）に反する。そのときζwn=ζe に変更して再度ρ5 を求める。ノッチフィルター効果（外乱除去比ζwn /ζew）を失うが、ρ5 が小さくなるのでωe5* を下げローパスフィルタ効果を得ることができる。すなわち（８５）式で
【０１９３】
【数１０４】
になり、Ｄc5Δb（ｓ）の３次推定器の特性に変更される。
【０１９４】
２．５ まとめ
２．４節の結果をまとめると、
・推定根のζe2*, ωe2* はζspec付近で、ρ2 にほぼ比例する。
・外乱モデルを含むρ5とωe5* との予測値は、外乱モデルを含まないρ2とωe2* によって求められる。
【０１９５】
７次特性根からρ5 を計算するとき、ρ5 とζe5*, ωe5* との初期値が適切であると求解が確実になる。
【０１９６】
以下、その求解算法について説明する。
【０１９７】
３．求解算法
推定係数ρを求める算法を示す。求められたρはパラメータ不確かさを印加した閉ループ系において、仕様Δspec, ζspecを満足する。算法は確実な求解法、簡単な処理と少ない計算量とのため、前章で導き出した性質を用いて、３つの段階により構築する。
【０１９８】
３．１ メインフロー
次の流れを採用する（図９）
（１）（５８）式でζe2*＝ζspec になるρ2，ωe2* を求める（Ｓ１０）。
（２）（８５）式でζe5*＝ζspec になるρ5，ωe5* を求めるために、（８６）式の予測値を初期値に設定する（Ｓ２０）。またωe5＞ωwn のとき（Ｓ３０でｙｅｓ）、ζwn＝ζe にして（Ｓ４０）、（２）を繰り返す。
（３）求めたρ5をρ（＝ωe／ωf）として制御系を構成する（Ｓ５０）。
【０１９９】
３．２ ステップ２ サブフロー１
ρ（ρ2、ρ5）はその性質を利用して２段階で求める（図１０参照）。
【０２００】
（１）ζspec を挟む範囲のρを求める：ζe* はρに対して単調増加になる性質を利用する。初期値は
【０２０１】
【数１０５】
とおく。ステップ３に初期値を送るとζe*，ωe* が戻される。反復ループ（添え字ｉ：回数）で
【０２０２】
【数１０６】
を更新する。上記にてｅζi＝ζei*−ζspec，ｅζi＞０のときブレークし、下限値と上限値とを設定する（図１０参照）。
【０２０３】
（２）次に、ζe*＝ζspec となるρを求める。ζe* はζspec 付近でρにほぼ直線近似になる性質を利用する。反復ループ（添え字ｉ：回数）で
【０２０４】
【数１０７】
を更新する。上記にて｜ｅζi｜＜ε（ε：微小項）のときブレークする。｜ｅζi｜≧εのときに、範囲を次式に置き換えて反複ループを実行する（図１０参照）。
【０２０５】
【数１０８】
【０２０６】
３．３ ステップ３
初期値を用いて次の処理を行い、ζe*とωe*とを取り出す。
（１）４次式は代数解法によって、最小の減衰係数に対応する２次式の係数を求める。
（２）７次式は任意の数値解法によって求める。具体的には、ベアストウ法 と組立除法とによって、２次式の係数を求め、残りの５次式を因数分解する場合は ニュートン−ラプソン法（初期値ωeo）と代数解法とを用いる。
【０２０７】
なお、上記の解法は公知のため説明を省く。
【０２０８】
４．計算結果
本発明を数値計算によって検証する。船体パラメータおよび制御パラメータを表３にまとめる。ただし設計パラメータについては表１を参照されたい。また、仕様は（４０）、（４２）式に従っている。
【０２０９】
【表３】
【０２１０】
４．１ ２次推定器
表３から以下のことが分かる。
・ρ2の変化は微分ゲインｆ2 のそれより小さいが、ωfＴsn に関係する。ρ2 は仕様から判断すると実用範囲の値である。
【０２１１】
４．２ ３次推定器
図１１を参照すると
・ａz /ｂ12：ωfＴsn ＝１ の場合は他の場合より１５〜３０％大きくなるが、ωfＴsn ＞１ の場合はほぼ近い値を示す。
・（−ａzｚ2）/ｂ02 ：近似値はａz /ｂ12 の９０％弱になり、ρo が大きくなると近似誤差が増加する。
・誤差ρ3~’，誤差ωe3~’：ρo≦０．２とすれば、両者とも５％以内の誤差に収まる。
【０２１２】
よって、３次推定器を２次推定器に近似し予測値を求める方法は妥当であり、予測値ρ3~’，ωe3~*’ は有効であると言える。
【０２１３】
４．３ ４次推定器
図１２（ωfＴsn ＝１．５８を利用）を参照すると
・ｂ14 /ｂ12 、ｋw：共にｒw ＝１で１に近く、ｒw ＝１０でほぼ１に収斂する。
・ζe4*’，ωe4*’：Ｂ4（ｓ）のみの場合（波浪成分Ｄew（ｓ），Ｄw（ｓ）を省く）は波浪成分を含む場合より変化量が小さい。
【０２１４】
図１３より、ρ4’ の予測値はζwn ＝０．１かつｒw が低減で誤差を生じるが、それ以外では誤差は小さいと見なせる。ωe4’ の予測値は一定のため、低い領域で誤差を生じる。
【０２１５】
これに対して、波浪成分による影響が支配的になり、予測値ρ4’ ，ωe4’ は実用精度をもつ。
【０２１６】
４．４ ５次推定器
・図１４は図１３にρo を追加した場合で、ほぼ同様の傾向を示し、オフセットは予測値とほぼ等価である。
・図１５はζwn＝０．１，ωfＴsn を変化させた場合で、ωfＴsn ＝１のｒw の低減で急増する。
【０２１７】
予測値ρ5’ ，ωe5’ はｒw の低減で波浪成分誤差をもつが、波浪パラメータ条件（破線より右側部分が条件を満足する領域）から予測値として実用精度をもつことが分かる。
【０２１８】
４．５ 計算時間の比較
本発明の求解算法の計算時間を標準算法と比較する。標準算法は既知の推定係数ρsol(
ｒw) から探索範囲［ρlow，ρhigh］＝ρsol (ｒw)［1/ｃp，ｃp］を与えて（Matlab のfminbnd を利用）、（１２）式の特性行列から特性根を求め、最小の減衰係数を取り出し、それが仕様と一致する推定係数を求める。計算結果を表４にまとめる。条件：図１４，ζwn＝０．１である。
【０２１９】
【表４】
【０２２０】
本発明の計算時間は０．５８２，０．５８１，０．５８５[ｓ]，平均０．５８３[ｓ] になり、標準算法の計算時間に比べ、１桁程度速いことが確認できた。
【０２２１】
５．まとめ
推定器を含んだ保針ループに対して、推定根に影響を与えるパラメータ不確かさを特定し、仕様を定め、推定根特性を仕様付近で把握することで、ρの予測値を求めた。外乱モデルなしの２次推定根特性を基に外乱モデルありの推定根特性を解析した。求解算法は数値解と推定根特性とを利用し予測値を初期値に採用した。
【０２２２】
数値計算によって、推定根特性と予測値との妥当性を確認し、求解算法の有効性を標準算法の計算時間と比べ確認した。
【符号の説明】
【０２２３】
１０ 自動操舵装置
１２ フィードバック制御器
１４ 操舵機
１６ 船体
１８ 制御対象
２２ 状態フィードバック制御器
２４ 推定器
【技術分野】
【０００１】
本発明は、推定器を有する船舶用自動操舵装置に関し、特にパラメータの不確かさの存在下で閉ループの安定性を図ることができる船舶用自動操舵装置に関する。
【背景技術】
【０００２】
船舶用自動操舵装置は、一定の設定方位に、ジャイロコンパスから検出される船首方位を追従させるために舵を制御する装置であり、保針と変針との機能を持つ。
【０００３】
その保針制御系は、大きな慣性力の船体を小さな制御力の操舵機によって制御し、かつ検出方位に含まれる外乱の大きさが方位と同等かそれ以上になる特徴をもつため、受動的な制御を実現するフィードバック制御器で制御される。図２に示したように、フィードバック制御器１２は、小さな比例ゲインを基にしたＰＤ制御を行う状態フィードバック制御器２２と、状態推定と外乱除去とを行う推定器２４から構成される。閉ループ安定性と外乱除去性への影響及び設計自由度は、実用上状態フィードバック制御器２２よりも推定器２４の方が大きい。
【０００４】
制御対象１８の船体モデルは、載荷の変化によるノミナル値のパラメータ不確かさを持ち、そのためノミナル値のモデルベースで構成する推定器２４は載荷の変化によって推定値に誤差を生じる。このため、推定値の状態フィードバックによる閉ループ安定性は、推定誤差によって劣化し、船体を蛇行航行させる（ヨーイング）おそれがある。
【０００５】
本発明者らは、かかる問題を解決するために、特許文献１において、船体パラメータのノミナル値のパラメータ不確かさを積極的に考慮に入れて、該パラメータ不確かさを起因とする推定誤差を小さくすることができるようにした船舶用自動操舵装置及びその推定器の設計方法を提案している。
ここでは、推定器の特性多項式として、
【０００６】
【数１】
と置き、λe2は、波浪モデルが仮に無いとしたときの船体モデルの状態量を推定するための特性多項式であり、ζe、ωeがそれぞれ船体モデルの状態量を推定するための減衰係数、固有周波数としたときに、ωeを操舵系周波数ωfのρ（＞１）倍(以下、ρ：推定係数)に設定し、ζeを１／√２に設定するようにして推定器を設計することにより、閉ループ安定性を低下させる要因となる推定角速度のパラメータ不確かさによる推定誤差を低減させることができる、ことを見出した。
【０００７】
さらに、特許文献１では、船体モデルと２つの外乱モデル（波浪モデル、舵角オフセットモデル）とからなる制御対象に対して、推定誤差が船体モデルのみのときの推定誤差と等価になるように推定器の固有周波数を外乱モデル特性から修正することを、提案している。これによれば、推定器の制御対象は船体モデルのみとなり、外乱モデルを無視することができるようになる。
【０００８】
また、特許文献２では、推定係数ρ＝ωe／ωfを、操舵ループの代表根に着目し、代表根の減衰係数とパラメータ不確かさとの関係から制御ゲインを設計している。即ち、ここでは推定器がない閉ループ系において、フィードバックゲインの設定、操舵系の特性およびパラメータ不確かさの影響を調べ、推定器がある閉ループ系において、特性多項式の導出、代表根の移動、代表根の減衰係数を調べて数値計算を実施した結果、パラメータ不確かさの影響を見積もることによって、推定器と操舵系との固有周波数を関係づける推定係数ρを設定している。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特許第４６０３８３２号公報
【特許文献２】特開２００７−１１８８２８号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかしながら、従来においては、パラメータ不確かさおよび外乱モデルの影響を間接的、近似的に扱うために、閉ループ安定性に多少の変動を持つという問題がある。即ち、推定ループにおいて閉ループ安定性を劣化させる場合があるが、従来の設計による装置では、推定ループにおけるパラメータ不確かさの影響を直接扱っていない。
【００１１】
また、船体モデルはパラメータ不確かさを複数もつために、特許文献２で採用したパラメータ不確かさ以外が船種、海象、載貨などにより無視できない状況を生じる場合があり、用いたパラメータ不確かさによる影響が減衰係数に間接的に作用するために局所的な安定性しか確保することができない場合がある、という問題がある。
【００１２】
本発明はかかる課題に鑑みなされたもので、外乱モデルを含んだ推定器と操舵系の閉ループにおいて、パラメータ不確かさに起因する推定誤差を低減することができる推定ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。
【００１３】
また、本発明の更なる目的は、推定ループの根を不安定にさせるパラメータ不確かさを特定し、それを用いて制御ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。
【課題を解決するための手段】
【００１４】
前述した目的を達成するために、請求項１記載の発明は、設定方位ψＲに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、
【００１５】
【数２】
前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、
【００１６】
【数３】
船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔを、以下で表したときに、
【００１７】
【数４】
【００１８】
パラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする。
【００１９】
請求項２記載の発明は、請求項１記載の船舶用自動操舵装置において、前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄc4Δを、
【００２０】
【数５】
としたときに、Ｄc4Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζe2＊とし、ζe2＊が第２仕様値（ζspec）の値を満足するωe2／ωf（＝ρ2）を求め、このρ2を用いた初期値とした数値解法によって前記ωe／ωf（＝ρ）を求めることを特徴とする。
【００２１】
請求項３記載の発明は、請求項２記載の船舶用自動操舵装置において、前記初期値は、
【００２２】
【数６】
で表されることを特徴とする。
【００２３】
請求項４記載の発明は、請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置において、前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさであることを特徴とする。
【００２４】
請求項５記載の発明は、請求項４記載の船舶用自動操舵装置において、前記特性多項式ＤcΔは、
【００２５】
【数７】
で表され、ここで、
Ｄw（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、
【００２６】
【数８】
と表され、Ｂ（ｓ）は、
【００２７】
【数９】
で表されることを特徴とする。
【発明の効果】
【００２８】
本発明によれば、外乱モデルとしての波浪モデル及び舵角オフセットモデルを含んだ推定器を含む閉ループの特性多項式に基づき、特性多項式に船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔからパラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することから、閉ループ安定性を増加させることができる。
【００２９】
また、前記パラメータ不確かさΔを、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさとすることで、閉ループ安定性に及ぼす影響の高いパラメータ不確かさに起因する推定根の不安定性を排除し、従来のものより信頼性を向上させることができる。
【００３０】
外乱モデルを含まない推定器の推定根から出発して、外乱モデルを含んだ推定器の推定根の予測値を求めて、特性多項式ＤcΔの特性根を求めることで、求解を高い収束性で短時間に行うことができるようになる。
【図面の簡単な説明】
【００３１】
【図１】本発明による船舶用自動操舵装置及び制御対象の全体構成を表すブロック図である。
【図１Ａ】船体モデルのブロック図である。ただし、舵感度時定数Ｔs3n項を省く。
【図２】閉ループ制御系を表すブロック図である。
【図３Ａ】２次推定器のブロック図である。
【図３Ｂ】３次推定器のブロック図である。
【図３Ｃ】４次推定器のブロック図である。
【図３Ｄ】５次推定器のブロック図である。
【図４】開ループ伝達関数による閉ループを表すブロック図である。
【図５】Ｄc4(s)の根軌跡を表すｓ平面である。
【図６】４次特性多項式の数値解から求めたρ2に対するζf2*、ωf2*、ζe2*、ωe2*とこれらの変化率の特性を表すグラフである。
【図７】仕様値付近での正規化したρ2に対するζe2*、ωe2*とこれらの変化率の特性を表すグラフである。
【図８】図６の結果から導いたρ2、ζe2*、ωe2*の正規化した場合の特性を表すグラフである。
【図９】本発明による求解法を表すフローチャートである。
【図１０】ρを求める原理を表す図である。
【図１１】３次推定器の推定根の特性を表す数値結果である。
【図１２】４次推定器の推定根の特性を表す数値結果である。
【図１３】４次推定器の推定根の特性を表す数値結果である。
【図１４】５次推定器の推定根の特性を表す数値結果である。
【図１５】５次推定器の推定根の特性を表す数値結果である。
【発明を実施するための形態】
【００３２】
以下、図面を用いて本発明の実施の形態を説明する。
【００３３】
１．１ 保針系の構成
本発明で対象とする船舶用自動操舵装置１０の保針モードは、図１に示すように操舵機１４及び船体１６を含む制御対象１８を設定方位ψRに船首方位ψを追従させるため、フィードバック制御器１２で偏差
【００３４】
【数１０】
から演算された舵角指令δcによって制御する。δcは操舵機１４の舵角δを経て船体に作用する。船体はパラメータ不確かさをもち外乱が印加する。操舵機１４は舵角および舵速度の振幅制限をもつが、保針時のδcおよび
の振幅が制限以下と見なせるので操舵機を省略する（δ＝δc）。なお設定方位ψRは保針時一定になるため、簡単化のためψR＝０とする。
【００３５】
１．２ 制御対象とパラメータ不確かさ
図１に示した制御対象１８は、操舵機１４を除き、船体モデルと外乱モデルとから構成される。
【００３６】
１．２．１ 船体モデル
船体モデルは、図１Ａに示される野本の線形モデル
【００３７】
【数１１】
を用いる。ここで、ｓはラプラス演算子を、ψは船首方位を、ｒを旋回角速度を、Ｋsn，Ｔsn，は操縦性指数（船体パラメータ）でそれぞれ旋回力ゲイン、時定数のノミナル値（添え字ｓ：船体、ｎ：ノミナル値）を示す。尚、（１）式の代わりに
【００３８】
【数１２】
を用いることもできるが、舵感度時定数Ｔs3n<Ｔsnであるため、以下、簡略化のために舵感度時定数Ｔs3nを０とする。
【００３９】
上式を状態空間表現で表すと
【００４０】
【数１３】
になる。ここで
【００４１】
（Ｔ：転置），ｒは旋回角速度（図１Ａ参照）を
【００４２】
【数１４】
【００４３】
【数１５】
Δa、Δbはそれぞれパラメータ不確かさで、実際値とノミナル値との差異によって
【００４４】
【数１６】
に定める。
【００４５】
１．２．２ 外乱モデル
外乱モデルは、図２に示すように、舵角オフセットモデルと波浪モデルとからなる。舵角オフセットモデルは風などに誘起された方位軸回りに作用するモーメントを舵角換算したものとする。波浪モデルは白色ノイズが入力した狭帯域フィルタの出力を方位換算したものとする。
【００４６】
波浪モデルを状態空間表現で表すと
【００４７】
【数１７】
【００４８】
【数１８】
になる。ここでδoは一定値とした舵角オフセット成分を、ｘw ＝ [ξ，ψw ]T ，ξ：変数，ψw は平均値ゼロの有色の波浪成分を、νは白色ノイズＮ（０，１）を
【００４９】
【数１９】
Ｋwn，ζwn，ωwn はそれぞれ波浪モデルのゲイン、減衰係数と中心周波数とを示す。波浪モデルはパラメータ不確かさを含まないとする。
【００５０】
１．３ 閉ループ制御系
１．３．１ フィードバック制御器
フィードバック制御器１２は、図２に示すように状態フィードバック制御器２２と推定器（オブザーバ）２４とからなる。閉ループ制御系は制御対象１８とフィードバック制御器１２とから構成し
【００５１】
【数２０】
になる。ここで
【００５２】
【数２１】
（^：推定値）を示し、
はフィードバックゲインを（ｆ1：比例ゲイン、ｆ2：微分ゲイン、１：舵角オフセットに対応），
は推定ゲインを
【００５３】
【数２２】
を示し、Ｏはゼロ行列を示す。推定器出力は
【００５４】
【数２３】
になる。
【００５５】
１．３．２ 閉ループ特性行列
閉ループ特性行列は、推定誤差を
【００５６】
【数２４】
を用いて上式より
【００５７】
【数２５】
【００５８】
【数２６】
になる。ここでＡcは特性行列を示す。船体モデルに状態フィードバックを行ったループをフィードバックループ（Ａc（１，１）、２次系）と呼び，制御対象に推定ゲインのフィードバックを行ったループを推定ループ（Ａc（２，２）、５次系）と呼ぶ。上式より要素Ａc（１，１）のフィードバックループによる固有値と、要素Ａc（２，２）の推定ループによる固有値は、ΔＡs，ΔＢsの影響を受ける。
【００５９】
１．４ ノミナル値による特性多項式
フィードバックループ、推定ループおよび波浪モデルのノミナル値の特性多項式を
【００６０】
【数２７】
【００６１】
【数２８】
【００６２】
【数２９】
とする。また、（１５）式は、
【００６３】
【数３０】
【００６４】
【数３１】
【００６５】
【数３２】
とする。ここでｄｅｔ：行列式、ｓ：ラプラス演算子、Ｉ：適当な単位行列、Ｄf：フィードバックループの特性多項式（添え字ｆ）、Ｄest：推定ループの特性多項式、Ｄw：波浪モデルの特性多項式（添え字ｗ）、Ｄe、Ｄew、Ｄeo：それぞれ船体モデル（添え字ｅ）、波浪モデル（添え字ｅｗ）と舵角オフセットモデル（添え字ｅｏ）とに対応する特性多項式、ζ、ω：２次系のそれぞれ減衰係数、固有周波数を示す。
【００６６】
１．５ 状態フィードバックゲイン
状態フィードバック制御器２２において、ノミナル値によるフィードバックループの状態フィードバックゲインを定める。状態フィードバックは図２に示すように、ＰＤ制御を用いており、
【００６７】
【数３３】
とおき、設計パラメータにｆ1とζfとを選ぶと、ｆ2とωfとはそれぞれ
【００６８】
【数３４】
【００６９】
【数３５】
として得られる。
【００７０】
１．６ 推定ゲイン
２次の推定ゲインは、図３Ａに示される外乱無しの船体モデルのみで設計される２次推定器におけるｋ1（＝ｋ11）、ｋ2（＝ｋ21）である。
【００７１】
３次の推定ゲインは、図３Ｂに示される舵角オフセットモデルを含めて設計される３次推定器におけるｋ1（＝ｋ13）、ｋ2（＝ｋ23）、ｋ5（＝ｋ53）である。
【００７２】
４次の推定ゲインは、図３Ｃに示される波浪モデルを含めて設計される４次推定器におけるｋ1（＝ｋ14）、ｋ2（＝ｋ24）、ｋ3（＝ｋ34）、ｋ4（＝ｋ４4）である。
【００７３】
５次の推定ゲインは、図３Ｄに示される波浪モデルと共に舵角オフセットモデルを含めて設計される５次推定器におけるｋ1（＝ｋ15）、ｋ2（＝ｋ25）、ｋ3（＝ｋ35）、ｋ4（＝ｋ45）、ｋ5（＝ｋ55）である。
【００７４】
５次推定器の特性多項式は、
【００７５】
【数３６】
になる。一方、設計する特性多項式は、（１５）式を展開して
【００７６】
【数３７】
になる。
よって、推定ゲインは上記２つの多項式の係数から次式より求まる。
【００７７】
【数３８】
ここで
【００７８】
【数３９】
を示す。
【００７９】
以降の設計に利用する推定ゲインを以下に示す。
【００８０】
【数４０】
【００８１】
【数４１】
【００８２】
【数４２】
【００８３】
また、舵感度時定数Ｔs3nを考慮した場合には、以下となる。
【００８４】
【数４３】
【００８５】
【数４４】
【００８６】
【数４５】
【００８７】
Ｄe（ｓ）のパラメータはフィードバックループのωfを基準に、
【００８８】
【数４６】
と定める。ここでρ＞１を推定係数と呼ぶ。ρを閉ループ安定性に基づいて以降で設計する。
【００８９】
１．７ 設計パラメータと制御パラメータ
特性多項式の設計パラメータと制御パラメータを、以下表１にまとめる。
【００９０】
【表１】
【００９１】
上記表より、Ｄe（ｓ）の減衰係数ζeと推定係数ρを与えれば、制御系が構成されることになる。尚、上記表から微分ゲインｆ2＞０、ζf＝１／√２と定めるため、（２３）式からωfＴsn＞１／√２になる。船体パラメータのノミナル値Ｋsn、Ｔsnは、船体の旋回の度に同定器によって同定され、更新される。波浪モデルの減衰係数ζwnと固有角周波数ωwnも周波数同定器によって同定される。
【００９２】
２ 設計
前章で検討したように、設計は、推定器の減衰係数ζeと推定係数ρに帰着する。まず、特性根で推定根を不安定にさせるパラメータ不確かさを特定し、仕様を定める。次いで、減衰係数ζeを推定ループ単体から推定誤差が低減する適正値を設定する。そして閉ループからρに対する仕様付近の推定根の特性を明らかにし、解の予測値を設定し、ρの求解算法を推定根特性及び予測値から算出する。
【００９３】
２．１ パラメータ不確かさ
２．１．１ 特性多項式の導出
（１２）式の特性行列の行列式である特性多項式は、ノミナル値による特性多項式とパラメータ不確かさによる多項式との和に整理される。特性多項式を因数分解すると特性根が得られる。すなわち
【００９４】
【数４７】
になる。推定器の次数ｉ（ｉ＝０、２、３、４、５）に対応した次数ｃｉの特性多項式は次式になる。
【００９５】
【数４８】
【００９６】
【数４９】
【００９７】
【数５０】
【００９８】
【数５１】
【００９９】
【数５２】
ここで、
【０１００】
【数５３】
を示す（ｉ：次数）。
【０１０１】
２．１．２ パラメータ不確かさの影響
パラメータ不確かさが閉ループ制御系に及ぼす影響を検討する。特性多項式Ｄc（ｓ）において、パラメータ不確かさΔa，Δbを開ループゲインと見なせば、根軌跡手法が適用でき、閉ループ安定性を検討できる。開ループ伝達関数ＧＨ（ｓ）を（３１）式より
【０１０２】
【数５４】
に定めると、閉ループ特性は図４に示すものと等価になる。根軌跡手法を用いて、推定器の次数と漸近線の本数（上式で分母の極数から分子のゼロ点数を差し引いた数）との関係を表２にまとめる。
【０１０３】
【表２】
【０１０４】
表２より２次から５次までの推定器の漸近線の本数と形状は同じであり、推定器に起因する閉ループ特性は、２次推定器を用いたＤc4（ｓ）を基本とすることができることが分かるから、推定器による閉ループ特性を把握するためには、Ｄc4（ｓ）を調べればよい。
【０１０５】
パラメータ不確かさΔa，Δb毎に対するＤc4（ｓ）の根軌跡をそれぞれ図５（ａ）及び（ｂ）に示す。同図で×は極を、○はゼロ点を示す。この導出の条件として、ωfＴsn ＝１．５８（Ｔsn ＝６２．５、Ｋsn ＝０．０４、ｆ1＝１）とする。
【０１０６】
図５から次のことが分かる。
・Δa＞０の場合はフィードバックループが不安定傾向になり、推定ループは間接的に影響する。よって、ζfの調整によって直接改善することができる。
・Δb＞０の場合は推定ループが不安定傾向になり、推定根が直接影響し、フィードバックループは間接的に影響する。
【０１０７】
従って、推定係数ρの設定はΔbを利用することがよいことが分かる。
【０１０８】
２．１．３ 設計仕様
推定係数ρが満足すべき仕様は、実際の制御対象の不確定要素、非線形性などを考慮して定める。
【０１０９】
具体的には、パラメータ不確かさの仕様Δspecは、Δbと対応するものとして、ある値（第１仕様値）に定める。そのため、パラメータ不確かさΔbの無次元量を導入する。
【０１１０】
【数５５】
ここで、（・）’は無次元量を意味する。
【０１１１】
そして、例えば、パラメータ不確かさの仕様Δspecは、パラメータ不確かさΔbがＫsn／Ｔsnの２倍に相当するものとして、
【０１１２】
【数５６】
に定める。
【０１１３】
また、閉ループ安定性の仕様ζspecは、特性根の中で最小の減衰係数に対応させて、ある値（第２仕様値）に決める。この値としては、例えば、０．４とすることができる。即ち、
【０１１４】
【数５７】
【０１１５】
【数５８】
である。ここで上付き文字(*) はパラメータ不確かさを含んだ場合を意味する。仕様ζspec＝０．４はステップ応答のオーバシュート２５%に相当する量で、負の実軸に対して偏角cos−１（０．４）＝６６．４［ｄｅｇ］になる（図５（ｂ）参照）。
【０１１６】
２．２ 減数係数ζe
外乱モデルを含まない２次推定器の減衰係数ζeはパラメータ不確かさに起因した推定誤差を低減させる値に選ぶ。
【０１１７】
船体モデルに対する推定ループは、（９）式及び図３Ａを参照すると、
【０１１８】
【数５９】
【０１１９】
【数６０】
になる。（４３）式から次式を得る。
【０１２０】
【数６１】
【０１２１】
上式及びそれを微分した式と実際値Ｋs、Ｔsの船体モデル（（１）式参照）
【０１２２】
【数６２】
とを（４４）式に代入すると、次式になる。
【０１２３】
【数６３】
ここで、Δe2、Δe1はパラメータ不確かさを表す。
【０１２４】
上式から伝達関数は、
【０１２５】
【数６４】
【０１２６】
【数６５】
になる。ここで、Ｄe（ｓ）は、（１７）式と等価である。
【０１２７】
また、（４８）式を（４５）式に代入すると、ｒの伝達関数は次式となる。
【０１２８】
【数６６】
【０１２９】
伝達関数Ｇψψ＾｜2（ｓ）、Ｇψｒ＾｜2（ｓ）は、Δe2、Δe1によって共に推定誤差を生じるが、閉ループ安定性により強く関与するＧψｒ＾｜2（ｓ）を用いて減衰係数ζeを定める。
【０１３０】
２．２．１ Δe1に関して
（５０）式と（２８）式とより、Δe1に関連する項を抽出すると、
【０１３１】
【数６７】
となる。Δe1を微小として、ζe*e1をΔe1に関する級数で近似すると、
【０１３２】
【数６８】
となる。これより、ζe*e1がΔe1による誤差を受けないようにするため、ｇ（ζe）＝０を解くと、次式を得る。
【０１３３】
【数６９】
【０１３４】
２．２．２ Δe2に関して
同様に、（５０）式と（２８）式とより、Δe2に関連する項を抽出すると、
【０１３５】
【数７０】
となる。Δe2を微小として、ζe*e2をΔe2に関する級数で近似すると、
【０１３６】
【数７１】
となる。これより、ζe*e2がΔe2による誤差を受けないようにするための条件は、
【０１３７】
【数７２】
となる。
【０１３８】
２．２．３ 減数係数ζeの決定
上記（５３）式、（５６）式に近い値を減数係数ζeの値として決定することが好ましい。具体的には、ωeＴsnが１より十分大きいものと仮定すると、（５３）式から減数係数ζeは、
【０１３９】
【数７３】
と定めることができ、この値は（５６）式からも、さほど離れない値であると期待できる。
【０１４０】
２．３ 推定係数ρ
推定係数ρは、（３１）式のパラメータ不確かさ項に依存する。まずは、船体パラメータに対応する（３３）式に基づいて、ρ2及び２次推定根の特性を把握する。そして、（３４）〜（３６）式のパラメータ不確かさ項の係数を、（３３）式のそれと比較することで、外乱モデルに起因するρｉ（ｉ＝３，４，５）及び３〜５次推定根の特性を把握する。
【０１４１】
手順としては、
・４次特性根とその２次推定根とを２次推定係数ρ2及びΔspecに関して解析し、ρ2と２次推定根との予測値を数値解法を利用して求める。
・５〜７次の特性多項式は、４次の特性多項式に基づいて解析し、５〜７次の特性根における推定係数と推定根との予測値をρ2と２次推定根とから求める。
【０１４２】
２．３．１ ２次推定器
（３１）式の４次特性多項式は、
【０１４３】
【数７４】
になる。ここで、Ｄ*（ｓ）は、パラメータ不確かさを含んだ多項式で、
【０１４４】
【数７５】
を示す。（３７）式、（２３）式及び（２７）式より、
【０１４５】
【数７６】
を示す。
【０１４６】
パラメータ不確かさは、（３９）式、（２２）式を用いて、
【０１４７】
【数７７】
になり、ΔbＢ2（ｓ）はｆ1が相殺され、ｆ1の影響を受けないことが分かる。
【０１４８】
次に、以降の３次推定器の５次特性多項式の解析の準備として、２次推定根の特性について検討することにする。２次推定根は、図５（ｂ）よりΔbの増加に伴って不安定側に移動し、その減衰係数ζe2*は４次特性根の中で最小になる。その特性を把握するために、ωfＴsn ＝１．５８（Ｔsn ＝６２．５、Ｋsn ＝０．０４、ｆ1＝１）において、Δb’を１，２（仕様値），３のそれぞれの場合に対して、ρ2を１〜１０まで増加させたときの数値解を求めて、ζf2*、ωf2*、ζe2*、ωe2*をそれぞれ求めると、図６に示すものとなる。この図６から、ρ2の増加に対応して、以下の傾向が把握される。
・ζe2*は増加し、ζe＝１／√２に漸近し、ωe2*は単調増加する。
・ζf2*、ωf2*は徐々に変化しなくなり、一定値に近づく。
・Δb’が大きくなるに従って、ζe2*、ωe2*の線形範囲が拡大する。
【０１４９】
よって、
を満足するρspec＝ρ2は存在し（図６の例では、ρ2＝４．５８）、そのときωspec＝ωe2*とおく。
【０１５０】
図７は、仕様付近の２次推定根の特性を示す。Δb’／Δspec＝１．０，１．１，１．２，１．３のそれぞれの場合の、ρ2’＝ρ2／ρspecを１〜１．４まで変化させたときの数値解を用いる。図７から以下のことが把握される。
・Δb’／Δspecが大きくなるほど、線形性は向上する。
・ζe2*，ωe2*はそれぞれのｄζe2*/ｄρ，ｄωe2*/ｄρがほぼ一定と扱え、ζspec付近でほぼ直線近似が成り立つ。
【０１５１】
【数７８】
【０１５２】
【数７９】
【０１５３】
図８は図７の結果から導いたρ2およびζe2*，ωe2*の正規化した特性を示す。横軸にΔb’ ／Δspecを、縦軸にζe2*’＝ζe2*／ζspec、ρ2’＝ρ2／ρspec、ωe2*’＝ωe2*／ωspecの無次元値を用いる。図８から以下のことが把握される。
・ζe2*’はΔb’ ／Δspecにほぼ反比例する。このことから、ζe2*’の予測値を
【０１５４】
【数８０】
と近似することができる。ここで、~：予測値を表す。
【０１５５】
近似誤差（＝（予測値）÷（計算値）−１）はΔb’／Δspec＝１．３で約２％になる。
・ρ2’、ωe2*’はほぼ同値でΔb’／Δspecにほぼ比例する。
【０１５６】
図８の読取から、仕様値に対するパラメータ不確かさであるΔb’ ／Δspec、に対応するρ2’、ωe2*’の予測値を次式で近似することができる。
【０１５７】
【数８１】
ここで近似誤差はΔb’／Δspec＝１．３のときρ~2’が約０．８％、ω~e2*’が約１．８％になる。
【０１５８】
このように、Δbの増加によるζe2*の劣化、ωe2*の変化およびρ2の増加は、Δb’／Δspecによる予測値から把握できることが分かる。
【０１５９】
２．３．２ ３次推定器
３次推定器は２次推定器に、舵角オフセット成分δoの推定を加えたもので、５次特性多項式は、（３４）式から、
【０１６０】
【数８２】
になる。ここで、Ｂ3（ｓ）は、（２８）式より、
【０１６１】
【数８３】
を示す。
【０１６２】
２次推定器との対比を行って、２次推定器からの変化の特性を把握することによって、３次推定器の解析を行う。まず、（６５）式のパラメータ不確かさの多項式を、
【０１６３】
【数８４】
の因数分解した形に置き換える。ここで、ζf ＝ ζe を用いて
【０１６４】
【数８５】
を示す。２次推定器で用いた具体的数値とρ2＝４．５８、ρo＝０．１から４ａzｃz／ｂz2＝０．３０５となり、４ａzｃz／ｂz2は小さいことが推定される。よって、
【０１６５】
【数８６】
なる近似を用いると（√（１＋ε）≒１＋ε／２ （０＜｜ε｜＜１）））、Ｚ（ｓ）の根として次式を得る（ｚ2＜ｚ1とする）。
【０１６６】
【数８７】
【０１６７】
【数８８】
【０１６８】
（６８）式において分母は１に近いため、ｚ1≒ωeoであり、しかも、Ｄeo3＊（ｓ）の特性根−ωeo3*は極−ωeoからゼロ点ｚ1の間にあるので、Ｄeo3*（ｓ）≒Ｄeo（ｓ）と近似することができる。これより（６５）式は次式になる。
【０１６９】
【数８９】
【０１７０】
【数９０】
【０１７１】
（７１）式を見ると、Ｄc43Δb（ｓ）は２次推定器の（５８）式のＤc4Δb（ｓ）と同じ形になっていることが分かる。そこで、両者のΔbの係数を比較するために、（５９）式のｓ1 係数ｂ12、ｓ0 係数ｂ02と（７１）式のｓ1 係数ａz、ｓ0 係数ａzｚ2とをそれぞれ比較すると、
【０１７２】
【数９１】
【０１７３】
【数９２】
【０１７４】
になる。（７２）式及び（７３）式の比がそれぞれ１より大きいことから、パラメータ不確かさΔbの係数は、３次推定器の５次特性多項式の方が、２次推定器の４次特性多項式に比べ、大きくなることが分かり、２次推定器におけるパラメータ不確かさΔbが等価的に大きくなったものとみなすことができる。そのため３次推定係数ρ3は、ρ2に対してパラメータ不確かさΔbの等価的増加分を修正する必要がある。Δbの係数比を具体的な数値を代入して計算すると（ρo＝０．１）
【０１７５】
【数９３】
になる。上式から大きい値のｓ1 係数を選ぶこととすると、２次推定器の４次特性多項式においてパラメータ不確かさがａz /ｂ12倍増加したことと等価になる。よって、前節の２次推定器の特性から（６４）式を用いて、ρ’3の予測値を換算すると
【０１７６】
【数９４】
と定めることができる。これよりρ3、ωe3*の予測値を次式に定めることができる。
【０１７７】
【数９５】
【０１７８】
２．３．３ ４次推定器
４次推定器は２次推定器に波浪成分ψwの推定を加えたもので、６次特性多項式は（３５）式から、
【０１７９】
【数９６】
になる。ここで、（２９）式より
【０１８０】
【数９７】
である。
【０１８１】
３次推定器での手法と同様に、２次推定器との対比を行って、４次推定器の解析を行う。そのため、Ｂ4の係数ｂ14と、Ｂ2の係数ｂ12との比較をすると、（７９）式から
【０１８２】
【数９８】
となる。詳細は省略するが、ｆ1Ｔsnωe22／ｂ12＞１であるため、上記式は、ｋwに関して右下がりの直線となる。
【０１８３】
波浪パラメータの中心角周波数比を
【０１８４】
【数９９】
と定め、ζwn＝０．１として数値計算すると
を得る。最悪ケースのｂ14 /ｂ12＝１．１２１は（７４）式のａz /ｂ12＝１．１３７とほぼ同等であるが、ｒw≫１でｂ14はｂ12に収束する。
【０１８５】
（７８）式でＤw(s)とＤew(s)とは次の条件のとき
【０１８６】
【数１００】
に近似できるので、２次推定器の特性と等価になり、４次推定係数ρ4≒ρ2，ωe4*≒ωe2*になる。上記の近似ができない場合、ρ4は減衰係数差（ζe−ζwn）に比例し、周波数比ｒwに反比例する傾向を持つとし、その予測値を次式に定める。
【０１８７】
【数１０１】
【０１８８】
２．３．４ ５次推定器
５次推定器は２次推定器に舵角オフセット成分δo、波浪成分ψwの推定を加えたもので、３次及び４次の推定根特性を踏襲するものとし、７次特性多項式は（３６）式から、次式となる。
【０１８９】
【数１０２】
【０１９０】
上式の予測値を３次、４次の推定根特性から次式に定める。
【０１９１】
【数１０３】
【０１９２】
２．４ 波浪パラメータ条件
Ｄc7Δb（ｓ）において、仕様を満たすがωe5 ＝ ρ5ωf＞ωwn になった場合、波浪パラメータ条件（表１）に反する。そのときζwn=ζe に変更して再度ρ5 を求める。ノッチフィルター効果（外乱除去比ζwn /ζew）を失うが、ρ5 が小さくなるのでωe5* を下げローパスフィルタ効果を得ることができる。すなわち（８５）式で
【０１９３】
【数１０４】
になり、Ｄc5Δb（ｓ）の３次推定器の特性に変更される。
【０１９４】
２．５ まとめ
２．４節の結果をまとめると、
・推定根のζe2*, ωe2* はζspec付近で、ρ2 にほぼ比例する。
・外乱モデルを含むρ5とωe5* との予測値は、外乱モデルを含まないρ2とωe2* によって求められる。
【０１９５】
７次特性根からρ5 を計算するとき、ρ5 とζe5*, ωe5* との初期値が適切であると求解が確実になる。
【０１９６】
以下、その求解算法について説明する。
【０１９７】
３．求解算法
推定係数ρを求める算法を示す。求められたρはパラメータ不確かさを印加した閉ループ系において、仕様Δspec, ζspecを満足する。算法は確実な求解法、簡単な処理と少ない計算量とのため、前章で導き出した性質を用いて、３つの段階により構築する。
【０１９８】
３．１ メインフロー
次の流れを採用する（図９）
（１）（５８）式でζe2*＝ζspec になるρ2，ωe2* を求める（Ｓ１０）。
（２）（８５）式でζe5*＝ζspec になるρ5，ωe5* を求めるために、（８６）式の予測値を初期値に設定する（Ｓ２０）。またωe5＞ωwn のとき（Ｓ３０でｙｅｓ）、ζwn＝ζe にして（Ｓ４０）、（２）を繰り返す。
（３）求めたρ5をρ（＝ωe／ωf）として制御系を構成する（Ｓ５０）。
【０１９９】
３．２ ステップ２ サブフロー１
ρ（ρ2、ρ5）はその性質を利用して２段階で求める（図１０参照）。
【０２００】
（１）ζspec を挟む範囲のρを求める：ζe* はρに対して単調増加になる性質を利用する。初期値は
【０２０１】
【数１０５】
とおく。ステップ３に初期値を送るとζe*，ωe* が戻される。反復ループ（添え字ｉ：回数）で
【０２０２】
【数１０６】
を更新する。上記にてｅζi＝ζei*−ζspec，ｅζi＞０のときブレークし、下限値と上限値とを設定する（図１０参照）。
【０２０３】
（２）次に、ζe*＝ζspec となるρを求める。ζe* はζspec 付近でρにほぼ直線近似になる性質を利用する。反復ループ（添え字ｉ：回数）で
【０２０４】
【数１０７】
を更新する。上記にて｜ｅζi｜＜ε（ε：微小項）のときブレークする。｜ｅζi｜≧εのときに、範囲を次式に置き換えて反複ループを実行する（図１０参照）。
【０２０５】
【数１０８】
【０２０６】
３．３ ステップ３
初期値を用いて次の処理を行い、ζe*とωe*とを取り出す。
（１）４次式は代数解法によって、最小の減衰係数に対応する２次式の係数を求める。
（２）７次式は任意の数値解法によって求める。具体的には、ベアストウ法 と組立除法とによって、２次式の係数を求め、残りの５次式を因数分解する場合は ニュートン−ラプソン法（初期値ωeo）と代数解法とを用いる。
【０２０７】
なお、上記の解法は公知のため説明を省く。
【０２０８】
４．計算結果
本発明を数値計算によって検証する。船体パラメータおよび制御パラメータを表３にまとめる。ただし設計パラメータについては表１を参照されたい。また、仕様は（４０）、（４２）式に従っている。
【０２０９】
【表３】
【０２１０】
４．１ ２次推定器
表３から以下のことが分かる。
・ρ2の変化は微分ゲインｆ2 のそれより小さいが、ωfＴsn に関係する。ρ2 は仕様から判断すると実用範囲の値である。
【０２１１】
４．２ ３次推定器
図１１を参照すると
・ａz /ｂ12：ωfＴsn ＝１ の場合は他の場合より１５〜３０％大きくなるが、ωfＴsn ＞１ の場合はほぼ近い値を示す。
・（−ａzｚ2）/ｂ02 ：近似値はａz /ｂ12 の９０％弱になり、ρo が大きくなると近似誤差が増加する。
・誤差ρ3~’，誤差ωe3~’：ρo≦０．２とすれば、両者とも５％以内の誤差に収まる。
【０２１２】
よって、３次推定器を２次推定器に近似し予測値を求める方法は妥当であり、予測値ρ3~’，ωe3~*’ は有効であると言える。
【０２１３】
４．３ ４次推定器
図１２（ωfＴsn ＝１．５８を利用）を参照すると
・ｂ14 /ｂ12 、ｋw：共にｒw ＝１で１に近く、ｒw ＝１０でほぼ１に収斂する。
・ζe4*’，ωe4*’：Ｂ4（ｓ）のみの場合（波浪成分Ｄew（ｓ），Ｄw（ｓ）を省く）は波浪成分を含む場合より変化量が小さい。
【０２１４】
図１３より、ρ4’ の予測値はζwn ＝０．１かつｒw が低減で誤差を生じるが、それ以外では誤差は小さいと見なせる。ωe4’ の予測値は一定のため、低い領域で誤差を生じる。
【０２１５】
これに対して、波浪成分による影響が支配的になり、予測値ρ4’ ，ωe4’ は実用精度をもつ。
【０２１６】
４．４ ５次推定器
・図１４は図１３にρo を追加した場合で、ほぼ同様の傾向を示し、オフセットは予測値とほぼ等価である。
・図１５はζwn＝０．１，ωfＴsn を変化させた場合で、ωfＴsn ＝１のｒw の低減で急増する。
【０２１７】
予測値ρ5’ ，ωe5’ はｒw の低減で波浪成分誤差をもつが、波浪パラメータ条件（破線より右側部分が条件を満足する領域）から予測値として実用精度をもつことが分かる。
【０２１８】
４．５ 計算時間の比較
本発明の求解算法の計算時間を標準算法と比較する。標準算法は既知の推定係数ρsol(
ｒw) から探索範囲［ρlow，ρhigh］＝ρsol (ｒw)［1/ｃp，ｃp］を与えて（Matlab のfminbnd を利用）、（１２）式の特性行列から特性根を求め、最小の減衰係数を取り出し、それが仕様と一致する推定係数を求める。計算結果を表４にまとめる。条件：図１４，ζwn＝０．１である。
【０２１９】
【表４】
【０２２０】
本発明の計算時間は０．５８２，０．５８１，０．５８５[ｓ]，平均０．５８３[ｓ] になり、標準算法の計算時間に比べ、１桁程度速いことが確認できた。
【０２２１】
５．まとめ
推定器を含んだ保針ループに対して、推定根に影響を与えるパラメータ不確かさを特定し、仕様を定め、推定根特性を仕様付近で把握することで、ρの予測値を求めた。外乱モデルなしの２次推定根特性を基に外乱モデルありの推定根特性を解析した。求解算法は数値解と推定根特性とを利用し予測値を初期値に採用した。
【０２２２】
数値計算によって、推定根特性と予測値との妥当性を確認し、求解算法の有効性を標準算法の計算時間と比べ確認した。
【符号の説明】
【０２２３】
１０ 自動操舵装置
１２ フィードバック制御器
１４ 操舵機
１６ 船体
１８ 制御対象
２２ 状態フィードバック制御器
２４ 推定器
【特許請求の範囲】
【請求項１】
設定方位ψＲに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、
【数１】
前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、
【数２】
船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔを、以下で表したときに、
【数３】
パラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする船舶用自動操舵装置。
【請求項２】
前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄc4Δを、
【数４】
としたときに、Ｄc4Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζe2＊とし、ζe2＊が第２仕様値（ζspec）の値を満足するωe2／ωf（＝ρ2）を求め、このρ2を用いた初期値とした数値解法によって前記ωe／ωf（＝ρ）を求めることを特徴とする請求項１記載の船舶用自動操舵装置。
【請求項３】
前記初期値は、
【数５】
で表されることを特徴とする請求項２記載の船舶用自動操舵装置。
【請求項４】
前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさであることを特徴とする請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置。
【請求項５】
前記特性多項式ＤcΔは、
【数６】
で表され、ここで、
Ｄw（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、
【数７】
と表され、Ｂ（ｓ）は、
【数８】
で表されることを特徴とする請求項４記載の船舶用自動操舵装置。
【請求項１】
設定方位ψＲに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、
【数１】
前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、
【数２】
船体パラメータのパラメータ不確かさΔを加味したときの特性多項式ＤcΔを、以下で表したときに、
【数３】
パラメータ不確かさΔが第１仕様値（Δspec）のときに、ＤcΔ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζe＊として、そのζe＊が第２仕様値（ζspec）の値を満足するωe／ωf（＝ρ）を求め、推定器の船体モデルの固有角周波数ωeを設定することを特徴とする船舶用自動操舵装置。
【請求項２】
前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄc4Δを、
【数４】
としたときに、Ｄc4Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζe2＊とし、ζe2＊が第２仕様値（ζspec）の値を満足するωe2／ωf（＝ρ2）を求め、このρ2を用いた初期値とした数値解法によって前記ωe／ωf（＝ρ）を求めることを特徴とする請求項１記載の船舶用自動操舵装置。
【請求項３】
前記初期値は、
【数５】
で表されることを特徴とする請求項２記載の船舶用自動操舵装置。
【請求項４】
前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫsnと時定数Ｔsnの比の不確かさであることを特徴とする請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置。
【請求項５】
前記特性多項式ＤcΔは、
【数６】
で表され、ここで、
Ｄw（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζwnと固有角周波数ωwnとで、
【数７】
と表され、Ｂ（ｓ）は、
【数８】
で表されることを特徴とする請求項４記載の船舶用自動操舵装置。
【図１】
【図１Ａ】
【図２】
【図３Ａ】
【図３Ｂ】
【図３Ｃ】
【図３Ｄ】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１Ａ】
【図２】
【図３Ａ】
【図３Ｂ】
【図３Ｃ】
【図３Ｄ】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【公開番号】特開２０１２−２１０８７５（Ｐ２０１２−２１０８７５Ａ）
【公開日】平成２４年１１月１日（２０１２．１１．１）
【国際特許分類】
【出願番号】特願２０１１−７７９７７（Ｐ２０１１−７７９７７）
【出願日】平成２３年３月３１日（２０１１．３．３１）
【出願人】（０００００３３８８）東京計器株式会社 (103)
【公開日】平成２４年１１月１日（２０１２．１１．１）
【国際特許分類】
【出願日】平成２３年３月３１日（２０１１．３．３１）
【出願人】（０００００３３８８）東京計器株式会社 (103)
[ Back to top ]