特許 5682009 船舶用自動操舵装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5682009

(24)【登録日】2015年1月23日

(45)【発行日】2015年3月11日

(54)【発明の名称】船舶用自動操舵装置

(51)【国際特許分類】

   B63H 25/04 20060101AFI20150219BHJP

【ＦＩ】

   B63H25/04 D

【請求項の数】5

【全頁数】36

(21)【出願番号】特願2011-77977(P2011-77977)

(22)【出願日】2011年3月31日

(65)【公開番号】特開2012-210875(P2012-210875A)

(43)【公開日】2012年11月1日

【審査請求日】2014年3月24日

(73)【特許権者】

【識別番号】000003388

【氏名又は名称】東京計器株式会社

(74)【代理人】

【識別番号】100101856

【弁理士】

【氏名又は名称】赤澤 日出夫

(72)【発明者】

【氏名】羽根 冬希

(72)【発明者】

【氏名】小山田 俊英

【審査官】 中村 泰二郎

(56)【参考文献】

【文献】 特開２００７−１１８８２８（ＪＰ，Ａ）

【文献】 特許第４６０３８３２（ＪＰ，Ｂ２）

【文献】 特開２００８−１３７５４５（ＪＰ，Ａ）

【文献】 米国特許出願公開第２００７／００５０１１３（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｂ６３Ｈ ２５／０４

Ｇ０８Ｇ ３／００

Ｇ０５Ｄ １／００

Ｇ０５Ｂ １３／０２

(57)【特許請求の範囲】

【請求項1】

設定方位ψ_Ｒに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、

【数1】

前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、

【数2】

船体パラメータのパラメータ不確かさΔを加味したときの特性多項式Ｄ_c^Δを、以下で表したときに、

【数3】

パラメータ不確かさΔが第１仕様値（Δ_spec）のときに、Ｄ_c^Δ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζ_e^＊として、そのζ_e^＊が第２仕様値（ζ_spec）の値を満足するω_e／ω_f（＝ρ）を求め、推定器の船体モデルの固有角周波数ω_eを設定することを特徴とする船舶用自動操舵装置。

【請求項2】

前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄ_c4^Δを、

【数4】

としたときに、Ｄ_c4^Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζ_e2^＊とし、ζ_e2^＊が第２仕様値（ζ_spec）の値を満足するω_e2／ω_f（＝ρ₂）を求め、このρ₂を用いた初期値とした数値解法によって前記ω_e／ω_f（＝ρ）を求めることを特徴とする請求項１記載の船舶用自動操舵装置。

【請求項3】

前記初期値は、

【数5】

で表されることを特徴とする請求項２記載の船舶用自動操舵装置。

【請求項4】

前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫ_snと時定数Ｔ_snの比の不確かさであることを特徴とする請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置。

【請求項5】

前記特性多項式Ｄ_c^Δは、

【数6】

で表され、ここで、
Ｄ_w（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζ_wnと固有角周波数ω_wnとで、

【数7】

と表され、Ｂ（ｓ）は、

【数8】

で表されることを特徴とする請求項４記載の船舶用自動操舵装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、推定器を有する船舶用自動操舵装置に関し、特にパラメータの不確かさの存在下で閉ループの安定性を図ることができる船舶用自動操舵装置に関する。

【背景技術】

【0002】

船舶用自動操舵装置は、一定の設定方位に、ジャイロコンパスから検出される船首方位を追従させるために舵を制御する装置であり、保針と変針との機能を持つ。

【0003】

その保針制御系は、大きな慣性力の船体を小さな制御力の操舵機によって制御し、かつ検出方位に含まれる外乱の大きさが方位と同等かそれ以上になる特徴をもつため、受動的な制御を実現するフィードバック制御器で制御される。図２に示したように、フィードバック制御器１２は、小さな比例ゲインを基にしたＰＤ制御を行う状態フィードバック制御器２２と、状態推定と外乱除去とを行う推定器２４から構成される。閉ループ安定性と外乱除去性への影響及び設計自由度は、実用上状態フィードバック制御器２２よりも推定器２４の方が大きい。

【0004】

制御対象１８の船体モデルは、載荷の変化によるノミナル値のパラメータ不確かさを持ち、そのためノミナル値のモデルベースで構成する推定器２４は載荷の変化によって推定値に誤差を生じる。このため、推定値の状態フィードバックによる閉ループ安定性は、推定誤差によって劣化し、船体を蛇行航行させる（ヨーイング）おそれがある。

【0005】

本発明者らは、かかる問題を解決するために、特許文献１において、船体パラメータのノミナル値のパラメータ不確かさを積極的に考慮に入れて、該パラメータ不確かさを起因とする推定誤差を小さくすることができるようにした船舶用自動操舵装置及びその推定器の設計方法を提案している。
ここでは、推定器の特性多項式として、

【0006】

【数1】

と置き、λ_e2は、波浪モデルが仮に無いとしたときの船体モデルの状態量を推定するための特性多項式であり、ζ_e、ω_eがそれぞれ船体モデルの状態量を推定するための減衰係数、固有周波数としたときに、ω_eを操舵系周波数ω_fのρ（＞１）倍(以下、ρ：推定係数)に設定し、ζ_eを１／√２に設定するようにして推定器を設計することにより、閉ループ安定性を低下させる要因となる推定角速度のパラメータ不確かさによる推定誤差を低減させることができる、ことを見出した。

【0007】

さらに、特許文献１では、船体モデルと２つの外乱モデル（波浪モデル、舵角オフセットモデル）とからなる制御対象に対して、推定誤差が船体モデルのみのときの推定誤差と等価になるように推定器の固有周波数を外乱モデル特性から修正することを、提案している。これによれば、推定器の制御対象は船体モデルのみとなり、外乱モデルを無視することができるようになる。

【0008】

また、特許文献２では、推定係数ρ＝ω_e／ω_fを、操舵ループの代表根に着目し、代表根の減衰係数とパラメータ不確かさとの関係から制御ゲインを設計している。即ち、ここでは推定器がない閉ループ系において、フィードバックゲインの設定、操舵系の特性およびパラメータ不確かさの影響を調べ、推定器がある閉ループ系において、特性多項式の導出、代表根の移動、代表根の減衰係数を調べて数値計算を実施した結果、パラメータ不確かさの影響を見積もることによって、推定器と操舵系との固有周波数を関係づける推定係数ρを設定している。

【先行技術文献】

【特許文献】

【0009】

【特許文献1】特許第４６０３８３２号公報

【特許文献2】特開２００７−１１８８２８号公報

【発明の概要】

【発明が解決しようとする課題】

【0010】

しかしながら、従来においては、パラメータ不確かさおよび外乱モデルの影響を間接的、近似的に扱うために、閉ループ安定性に多少の変動を持つという問題がある。即ち、推定ループにおいて閉ループ安定性を劣化させる場合があるが、従来の設計による装置では、推定ループにおけるパラメータ不確かさの影響を直接扱っていない。

【0011】

また、船体モデルはパラメータ不確かさを複数もつために、特許文献２で採用したパラメータ不確かさ以外が船種、海象、載貨などにより無視できない状況を生じる場合があり、用いたパラメータ不確かさによる影響が減衰係数に間接的に作用するために局所的な安定性しか確保することができない場合がある、という問題がある。

【0012】

本発明はかかる課題に鑑みなされたもので、外乱モデルを含んだ推定器と操舵系の閉ループにおいて、パラメータ不確かさに起因する推定誤差を低減することができる推定ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。

【0013】

また、本発明の更なる目的は、推定ループの根を不安定にさせるパラメータ不確かさを特定し、それを用いて制御ゲインを設定した推定器を有する船舶用自動操舵装置を提供することをその目的とする。

【課題を解決するための手段】

【0014】

前述した目的を達成するために、請求項１記載の発明は、設定方位ψ_Ｒに検出方位を追従させるべく制御対象を制御する船舶用自動操舵装置であって、制御対象からの検出方位から外乱を除去した推定値を出力する推定器と、推定器からの推定値と設定方位から得られる値に対してフィードバックゲインを乗じて命令舵角を出力する状態フィードバック制御器とを有するフィードバック制御器を備え、命令舵角に応じて舵角を経て船体に作用させる操舵機及び船体を含む制御対象と共に閉ループを構成する船舶用自動操舵装置において、
前記推定器の特性多項式は、船体モデルと波浪モデルと舵角オフセットモデルとを推定する多項式からなるものとし、それぞれ以下で表し、

【0015】

【数2】

前記閉ループの特性多項式は、フィードバックループの特性多項式と、推定器の特性多項式から以下で表し、

【0016】

【数3】

船体パラメータのパラメータ不確かさΔを加味したときの特性多項式Ｄ_c^Δを、以下で表したときに、

【0017】

【数4】

【0018】

パラメータ不確かさΔが第１仕様値（Δ_spec）のときに、Ｄ_c^Δ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζ_e^＊として、そのζ_e^＊が第２仕様値（ζ_spec）の値を満足するω_e／ω_f（＝ρ）を求め、推定器の船体モデルの固有角周波数ω_eを設定することを特徴とする。

【0019】

請求項２記載の発明は、請求項１記載の船舶用自動操舵装置において、前記推定器を船体モデルのみとしたときの前記閉ループの特性多項式Ｄ_c4^Δを、

【0020】

【数5】

としたときに、Ｄ_c4^Δ（ｓ）の特性根のうちの最も小さい減衰係数ζをζ_e2^＊とし、ζ_e2^＊が第２仕様値（ζ_spec）の値を満足するω_e2／ω_f（＝ρ₂）を求め、このρ₂を用いた初期値とした数値解法によって前記ω_e／ω_f（＝ρ）を求めることを特徴とする。

【0021】

請求項３記載の発明は、請求項２記載の船舶用自動操舵装置において、前記初期値は、

【0022】

【数6】

で表されることを特徴とする。

【0023】

請求項４記載の発明は、請求項１ないし３のいずれか１項に記載の船舶用自動操舵装置において、前記パラメータ不確かさΔは、船体パラメータの旋回力ゲインＫ_snと時定数Ｔ_snの比の不確かさであることを特徴とする。

【0024】

請求項５記載の発明は、請求項４記載の船舶用自動操舵装置において、前記特性多項式Ｄ_c^Δは、

【0025】

【数7】

で表され、ここで、
Ｄ_w（ｓ）は、波浪モデルの特性多項式であり、減衰係数ζ_wnと固有角周波数ω_wnとで、

【0026】

【数8】

と表され、Ｂ（ｓ）は、

【0027】

【数9】

で表されることを特徴とする。

【発明の効果】

【0028】

本発明によれば、外乱モデルとしての波浪モデル及び舵角オフセットモデルを含んだ推定器を含む閉ループの特性多項式に基づき、特性多項式に船体パラメータのパラメータ不確かさΔを加味したときの特性多項式Ｄ_c^Δからパラメータ不確かさΔが第１仕様値（Δ_spec）のときに、Ｄ_c^Δ（ｓ）の特性根から得られる減衰係数のうちの最も小さい減衰係数ζをζ_e^＊として、そのζ_e^＊が第２仕様値（ζ_spec）の値を満足するω_e／ω_f（＝ρ）を求め、推定器の船体モデルの固有角周波数ω_eを設定することから、閉ループ安定性を増加させることができる。

【0029】

また、前記パラメータ不確かさΔを、船体パラメータの旋回力ゲインＫ_snと時定数Ｔ_snの比の不確かさとすることで、閉ループ安定性に及ぼす影響の高いパラメータ不確かさに起因する推定根の不安定性を排除し、従来のものより信頼性を向上させることができる。

【0030】

外乱モデルを含まない推定器の推定根から出発して、外乱モデルを含んだ推定器の推定根の予測値を求めて、特性多項式Ｄ_c^Δの特性根を求めることで、求解を高い収束性で短時間に行うことができるようになる。

【図面の簡単な説明】

【0031】

【図1】本発明による船舶用自動操舵装置及び制御対象の全体構成を表すブロック図である。

【図1A】船体モデルのブロック図である。ただし、舵感度時定数Ｔ_s3n項を省く。

【図2】閉ループ制御系を表すブロック図である。

【図3A】２次推定器のブロック図である。

【図3B】３次推定器のブロック図である。

【図3C】４次推定器のブロック図である。

【図3D】５次推定器のブロック図である。

【図4】開ループ伝達関数による閉ループを表すブロック図である。

【図5】Ｄ_c4(s)の根軌跡を表すｓ平面である。

【図6】４次特性多項式の数値解から求めたρ₂に対するζ_f2^*、ω_f2^*_、ζ_e2^*、ω_e2^*とこれらの変化率の特性を表すグラフである。

【図7】仕様値付近での正規化したρ₂に対するζ_e2^*、ω_e2^*とこれらの変化率の特性を表すグラフである。

【図8】図６の結果から導いたρ₂、ζ_e2^*、ω_e2^*の正規化した場合の特性を表すグラフである。

【図9】本発明による求解法を表すフローチャートである。

【図10】ρを求める原理を表す図である。

【図11】３次推定器の推定根の特性を表す数値結果である。

【図12】４次推定器の推定根の特性を表す数値結果である。

【図13】４次推定器の推定根の特性を表す数値結果である。

【図14】５次推定器の推定根の特性を表す数値結果である。

【図15】５次推定器の推定根の特性を表す数値結果である。

【発明を実施するための形態】

【0032】

以下、図面を用いて本発明の実施の形態を説明する。

【0033】

１．１ 保針系の構成
本発明で対象とする船舶用自動操舵装置１０の保針モードは、図１に示すように操舵機１４及び船体１６を含む制御対象１８を設定方位ψ_Rに船首方位ψを追従させるため、フィードバック制御器１２で偏差

【0034】

【数10】

から演算された舵角指令δ_cによって制御する。δ_cは操舵機１４の舵角δを経て船体に作用する。船体はパラメータ不確かさをもち外乱が印加する。操舵機１４は舵角および舵速度の振幅制限をもつが、保針時のδ_cおよび

の振幅が制限以下と見なせるので操舵機を省略する（δ＝δ_c）。なお設定方位ψ_Rは保針時一定になるため、簡単化のためψ_R＝０とする。

【0035】

１．２ 制御対象とパラメータ不確かさ
図１に示した制御対象１８は、操舵機１４を除き、船体モデルと外乱モデルとから構成される。

【0036】

１．２．１ 船体モデル
船体モデルは、図１Ａに示される野本の線形モデル

【0037】

【数11】

を用いる。ここで、ｓはラプラス演算子を、ψは船首方位を、ｒを旋回角速度を、Ｋ_sn，Ｔ_sn，は操縦性指数（船体パラメータ）でそれぞれ旋回力ゲイン、時定数のノミナル値（添え字ｓ：船体、ｎ：ノミナル値）を示す。尚、（１）式の代わりに

【0038】

【数12】

を用いることもできるが、舵感度時定数Ｔ_s3n<Ｔ_snであるため、以下、簡略化のために舵感度時定数Ｔ_s3nを０とする。

【0039】

上式を状態空間表現で表すと

【0040】

【数13】

になる。ここで

【0041】

（Ｔ：転置），ｒは旋回角速度（図１Ａ参照）を

【0042】

【数14】

【0043】

【数15】

Δ_a、Δ_bはそれぞれパラメータ不確かさで、実際値とノミナル値との差異によって

【0044】

【数16】

に定める。

【0045】

１．２．２ 外乱モデル
外乱モデルは、図２に示すように、舵角オフセットモデルと波浪モデルとからなる。舵角オフセットモデルは風などに誘起された方位軸回りに作用するモーメントを舵角換算したものとする。波浪モデルは白色ノイズが入力した狭帯域フィルタの出力を方位換算したものとする。

【0046】

波浪モデルを状態空間表現で表すと

【0047】

【数17】

【0048】

【数18】

になる。ここでδ_oは一定値とした舵角オフセット成分を、ｘ_w ＝ [ξ，ψ_w ]^T ，ξ：変数，ψ_w は平均値ゼロの有色の波浪成分を、νは白色ノイズＮ（０，１）を

【0049】

【数19】

Ｋ_wn，ζ_wn，ω_wn はそれぞれ波浪モデルのゲイン、減衰係数と中心周波数とを示す。波浪モデルはパラメータ不確かさを含まないとする。

【0050】

１．３ 閉ループ制御系
１．３．１ フィードバック制御器
フィードバック制御器１２は、図２に示すように状態フィードバック制御器２２と推定器（オブザーバ）２４とからなる。閉ループ制御系は制御対象１８とフィードバック制御器１２とから構成し

【0051】

【数20】

になる。ここで

【0052】

【数21】

（^：推定値）を示し、

はフィードバックゲインを（ｆ₁：比例ゲイン、ｆ₂：微分ゲイン、１：舵角オフセットに対応），

は推定ゲインを

【0053】

【数22】

を示し、Ｏはゼロ行列を示す。推定器出力は

【0054】

【数23】

になる。

【0055】

１．３．２ 閉ループ特性行列
閉ループ特性行列は、推定誤差を

【0056】

【数24】

を用いて上式より

【0057】

【数25】

【0058】

【数26】

になる。ここでＡ_cは特性行列を示す。船体モデルに状態フィードバックを行ったループをフィードバックループ（Ａ_c（１，１）、２次系）と呼び，制御対象に推定ゲインのフィードバックを行ったループを推定ループ（Ａ_c（２，２）、５次系）と呼ぶ。上式より要素Ａ_c（１，１）のフィードバックループによる固有値と、要素Ａ_c（２，２）の推定ループによる固有値は、ΔＡ_s，ΔＢ_sの影響を受ける。

【0059】

１．４ ノミナル値による特性多項式
フィードバックループ、推定ループおよび波浪モデルのノミナル値の特性多項式を

【0060】

【数27】

【0061】

【数28】

【0062】

【数29】

とする。また、（１５）式は、

【0063】

【数30】

【0064】

【数31】

【0065】

【数32】

とする。ここでｄｅｔ：行列式、ｓ：ラプラス演算子、Ｉ：適当な単位行列、Ｄ_f：フィードバックループの特性多項式（添え字ｆ）、Ｄ_est：推定ループの特性多項式、Ｄ_w：波浪モデルの特性多項式（添え字ｗ）、Ｄ_e、Ｄ_ew、Ｄ_eo：それぞれ船体モデル（添え字ｅ）、波浪モデル（添え字ｅｗ）と舵角オフセットモデル（添え字ｅｏ）とに対応する特性多項式、ζ、ω：２次系のそれぞれ減衰係数、固有周波数を示す。

【0066】

１．５ 状態フィードバックゲイン
状態フィードバック制御器２２において、ノミナル値によるフィードバックループの状態フィードバックゲインを定める。状態フィードバックは図２に示すように、ＰＤ制御を用いており、

【0067】

【数33】

とおき、設計パラメータにｆ₁とζ_fとを選ぶと、ｆ₂とω_fとはそれぞれ

【0068】

【数34】

【0069】

【数35】

として得られる。

【0070】

１．６ 推定ゲイン
２次の推定ゲインは、図３Ａに示される外乱無しの船体モデルのみで設計される２次推定器におけるｋ₁（＝ｋ₁₁）、ｋ₂（＝ｋ₂₁）である。

【0071】

３次の推定ゲインは、図３Ｂに示される舵角オフセットモデルを含めて設計される３次推定器におけるｋ₁（＝ｋ₁₃）、ｋ₂（＝ｋ₂₃）、ｋ₅（＝ｋ₅₃）である。

【0072】

４次の推定ゲインは、図３Ｃに示される波浪モデルを含めて設計される４次推定器におけるｋ₁（＝ｋ₁₄）、ｋ₂（＝ｋ₂₄）、ｋ₃（＝ｋ₃₄）、ｋ₄（＝ｋ_４4）である。

【0073】

５次の推定ゲインは、図３Ｄに示される波浪モデルと共に舵角オフセットモデルを含めて設計される５次推定器におけるｋ₁（＝ｋ₁₅）、ｋ₂（＝ｋ₂₅）、ｋ₃（＝ｋ₃₅）、ｋ₄（＝ｋ₄₅）、ｋ₅（＝ｋ₅₅）である。

【0074】

５次推定器の特性多項式は、

【0075】

【数36】

になる。一方、設計する特性多項式は、（１５）式を展開して

【0076】

【数37】

になる。
よって、推定ゲインは上記２つの多項式の係数から次式より求まる。

【0077】

【数38】

ここで

【0078】

【数39】

を示す。

【0079】

以降の設計に利用する推定ゲインを以下に示す。

【0080】

【数40】

【0081】

【数41】

【0082】

【数42】

【0083】

また、舵感度時定数Ｔ_s3nを考慮した場合には、以下となる。

【0084】

【数43】

【0085】

【数44】

【0086】

【数45】

【0087】

Ｄ_e（ｓ）のパラメータはフィードバックループのω_fを基準に、

【0088】

【数46】

と定める。ここでρ＞１を推定係数と呼ぶ。ρを閉ループ安定性に基づいて以降で設計する。

【0089】

１．７ 設計パラメータと制御パラメータ
特性多項式の設計パラメータと制御パラメータを、以下表１にまとめる。

【0090】

【表1】

【0091】

上記表より、Ｄ_e（ｓ）の減衰係数ζ_eと推定係数ρを与えれば、制御系が構成されることになる。尚、上記表から微分ゲインｆ₂＞０、ζ_f＝１／√２と定めるため、（２３）式からω_fＴ_sn＞１／√２になる。船体パラメータのノミナル値Ｋ_sn、Ｔ_snは、船体の旋回の度に同定器によって同定され、更新される。波浪モデルの減衰係数ζ_wnと固有角周波数ω_wnも周波数同定器によって同定される。

【0092】

２ 設計
前章で検討したように、設計は、推定器の減衰係数ζ_eと推定係数ρに帰着する。まず、特性根で推定根を不安定にさせるパラメータ不確かさを特定し、仕様を定める。次いで、減衰係数ζ_eを推定ループ単体から推定誤差が低減する適正値を設定する。そして閉ループからρに対する仕様付近の推定根の特性を明らかにし、解の予測値を設定し、ρの求解算法を推定根特性及び予測値から算出する。

【0093】

２．１ パラメータ不確かさ
２．１．１ 特性多項式の導出
（１２）式の特性行列の行列式である特性多項式は、ノミナル値による特性多項式とパラメータ不確かさによる多項式との和に整理される。特性多項式を因数分解すると特性根が得られる。すなわち

【0094】

【数47】

になる。推定器の次数ｉ（ｉ＝０、２、３、４、５）に対応した次数ｃｉの特性多項式は次式になる。

【0095】

【数48】

【0096】

【数49】

【0097】

【数50】

【0098】

【数51】

【0099】

【数52】

ここで、

【0100】

【数53】

を示す（ｉ：次数）。

【0101】

２．１．２ パラメータ不確かさの影響
パラメータ不確かさが閉ループ制御系に及ぼす影響を検討する。特性多項式Ｄ_c（ｓ）において、パラメータ不確かさΔ_a，Δ_bを開ループゲインと見なせば、根軌跡手法が適用でき、閉ループ安定性を検討できる。開ループ伝達関数ＧＨ（ｓ）を（３１）式より

【0102】

【数54】

に定めると、閉ループ特性は図４に示すものと等価になる。根軌跡手法を用いて、推定器の次数と漸近線の本数（上式で分母の極数から分子のゼロ点数を差し引いた数）との関係を表２にまとめる。

【0103】

【表2】

【0104】

表２より２次から５次までの推定器の漸近線の本数と形状は同じであり、推定器に起因する閉ループ特性は、２次推定器を用いたＤ_c4（ｓ）を基本とすることができることが分かるから、推定器による閉ループ特性を把握するためには、Ｄ_c4（ｓ）を調べればよい。

【0105】

パラメータ不確かさΔ_a，Δ_b毎に対するＤ_c4（ｓ）の根軌跡をそれぞれ図５（ａ）及び（ｂ）に示す。同図で×は極を、○はゼロ点を示す。この導出の条件として、ω_fＴ_sn ＝１．５８（Ｔ_sn ＝６２．５、Ｋ_sn ＝０．０４、ｆ₁＝１）とする。

【0106】

図５から次のことが分かる。
・Δ_a＞０の場合はフィードバックループが不安定傾向になり、推定ループは間接的に影響する。よって、ζ_fの調整によって直接改善することができる。
・Δ_b＞０の場合は推定ループが不安定傾向になり、推定根が直接影響し、フィードバックループは間接的に影響する。

【0107】

従って、推定係数ρの設定はΔ_bを利用することがよいことが分かる。

【0108】

２．１．３ 設計仕様
推定係数ρが満足すべき仕様は、実際の制御対象の不確定要素、非線形性などを考慮して定める。

【0109】

具体的には、パラメータ不確かさの仕様Δ_specは、Δ_bと対応するものとして、ある値（第１仕様値）に定める。そのため、パラメータ不確かさΔ_bの無次元量を導入する。

【0110】

【数55】

ここで、（・）’は無次元量を意味する。

【0111】

そして、例えば、パラメータ不確かさの仕様Δ_specは、パラメータ不確かさΔ_bがＫ_sn／Ｔ_snの２倍に相当するものとして、

【0112】

【数56】

に定める。

【0113】

また、閉ループ安定性の仕様ζ_specは、特性根の中で最小の減衰係数に対応させて、ある値（第２仕様値）に決める。この値としては、例えば、０．４とすることができる。即ち、

【0114】

【数57】

【0115】

【数58】

である。ここで上付き文字(*) はパラメータ不確かさを含んだ場合を意味する。仕様ζ_spec＝０．４はステップ応答のオーバシュート２５%に相当する量で、負の実軸に対して偏角cos^−１（０．４）＝６６．４［ｄｅｇ］になる（図５（ｂ）参照）。

【0116】

２．２ 減数係数ζ_e
外乱モデルを含まない２次推定器の減衰係数ζ_eはパラメータ不確かさに起因した推定誤差を低減させる値に選ぶ。

【0117】

船体モデルに対する推定ループは、（９）式及び図３Ａを参照すると、

【0118】

【数59】

【0119】

【数60】

になる。（４３）式から次式を得る。

【0120】

【数61】

【0121】

上式及びそれを微分した式と実際値Ｋ_s、Ｔ_sの船体モデル（（１）式参照）

【0122】

【数62】

とを（４４）式に代入すると、次式になる。

【0123】

【数63】

ここで、Δ_e2、Δ_e1はパラメータ不確かさを表す。

【0124】

上式から伝達関数は、

【0125】

【数64】

【0126】

【数65】

になる。ここで、Ｄ_e（ｓ）は、（１７）式と等価である。

【0127】

また、（４８）式を（４５）式に代入すると、ｒの伝達関数は次式となる。

【0128】

【数66】

【0129】

伝達関数Ｇ_ψ^ψ＾｜₂（ｓ）、Ｇ_ψ^ｒ＾｜₂（ｓ）は、Δ_e2、Δ_e1によって共に推定誤差を生じるが、閉ループ安定性により強く関与するＧ_ψ^ｒ＾｜₂（ｓ）を用いて減衰係数ζ_eを定める。

【0130】

２．２．１ Δ_e1に関して
（５０）式と（２８）式とより、Δ_e1に関連する項を抽出すると、

【0131】

【数67】

となる。Δ_e1を微小として、ζ_e^*_e1をΔ_e1に関する級数で近似すると、

【0132】

【数68】

となる。これより、ζ_e^*_e1がΔ_e1による誤差を受けないようにするため、ｇ（ζ_e）＝０を解くと、次式を得る。

【0133】

【数69】

【0134】

２．２．２ Δ_e2に関して
同様に、（５０）式と（２８）式とより、Δ_e2に関連する項を抽出すると、

【0135】

【数70】

となる。Δ_e2を微小として、ζ_e^*_e2をΔ_e2に関する級数で近似すると、

【0136】

【数71】

となる。これより、ζ_e^*_e2がΔ_e2による誤差を受けないようにするための条件は、

【0137】

【数72】

となる。

【0138】

２．２．３ 減数係数ζ_eの決定
上記（５３）式、（５６）式に近い値を減数係数ζ_eの値として決定することが好ましい。具体的には、ω_eＴ_snが１より十分大きいものと仮定すると、（５３）式から減数係数ζ_eは、

【0139】

【数73】

と定めることができ、この値は（５６）式からも、さほど離れない値であると期待できる。

【0140】

２．３ 推定係数ρ
推定係数ρは、（３１）式のパラメータ不確かさ項に依存する。まずは、船体パラメータに対応する（３３）式に基づいて、ρ₂及び２次推定根の特性を把握する。そして、（３４）〜（３６）式のパラメータ不確かさ項の係数を、（３３）式のそれと比較することで、外乱モデルに起因するρ_ｉ（ｉ＝３，４，５）及び３〜５次推定根の特性を把握する。

【0141】

手順としては、
・４次特性根とその２次推定根とを２次推定係数ρ₂及びΔ_specに関して解析し、ρ₂と２次推定根との予測値を数値解法を利用して求める。
・５〜７次の特性多項式は、４次の特性多項式に基づいて解析し、５〜７次の特性根における推定係数と推定根との予測値をρ₂と２次推定根とから求める。

【0142】

２．３．１ ２次推定器
（３１）式の４次特性多項式は、

【0143】

【数74】

になる。ここで、Ｄ^*（ｓ）は、パラメータ不確かさを含んだ多項式で、

【0144】

【数75】

を示す。（３７）式、（２３）式及び（２７）式より、

【0145】

【数76】

を示す。

【0146】

パラメータ不確かさは、（３９）式、（２２）式を用いて、

【0147】

【数77】

になり、Δ_bＢ₂（ｓ）はｆ₁が相殺され、ｆ₁の影響を受けないことが分かる。

【0148】

次に、以降の３次推定器の５次特性多項式の解析の準備として、２次推定根の特性について検討することにする。２次推定根は、図５（ｂ）よりΔ_bの増加に伴って不安定側に移動し、その減衰係数ζ_e2^*は４次特性根の中で最小になる。その特性を把握するために、ω_fＴ_sn ＝１．５８（Ｔ_sn ＝６２．５、Ｋ_sn ＝０．０４、ｆ₁＝１）において、Δ_b’を１，２（仕様値），３のそれぞれの場合に対して、ρ₂を１〜１０まで増加させたときの数値解を求めて、ζ_f2^*、ω_f2^*_、ζ_e2^*、ω_e2^*をそれぞれ求めると、図６に示すものとなる。この図６から、ρ₂の増加に対応して、以下の傾向が把握される。
・ζ_e2^*は増加し、ζ_e＝１／√２に漸近し、ω_e2^*は単調増加する。
・ζ_f2^*、ω_f2^*は徐々に変化しなくなり、一定値に近づく。
・Δ_b’が大きくなるに従って、ζ_e2^*、ω_e2^*の線形範囲が拡大する。

【0149】

よって、

を満足するρ_spec＝ρ₂は存在し（図６の例では、ρ₂＝４．５８）、そのときω_spec＝ω_e2^*とおく。

【0150】

図７は、仕様付近の２次推定根の特性を示す。Δ_b’／Δ_spec＝１．０，１．１，１．２，１．３のそれぞれの場合の、ρ₂’＝ρ₂／ρ_specを１〜１．４まで変化させたときの数値解を用いる。図７から以下のことが把握される。
・Δ_b’／Δ_specが大きくなるほど、線形性は向上する。
・ζ_e2^*，ω_e2^*はそれぞれのｄζ_e2^*/ｄρ，ｄω_e2^*/ｄρがほぼ一定と扱え、ζ_spec付近でほぼ直線近似が成り立つ。

【0151】

【数78】

【0152】

【数79】

【0153】

図８は図７の結果から導いたρ₂およびζ_e2^*，ω_e2^*の正規化した特性を示す。横軸にΔ_b^’ ／Δ_specを、縦軸にζ_e2^*’＝ζ_e2^*／ζ_spec、ρ₂’＝ρ₂／ρ_spec、ω_e2^*’＝ω_e2^*／ω_specの無次元値を用いる。図８から以下のことが把握される。
・ζ_e2^*’はΔ_b’ ／Δ_specにほぼ反比例する。このことから、ζ_e2^*’の予測値を

【0154】

【数80】

と近似することができる。ここで、~：予測値を表す。

【0155】

近似誤差（＝（予測値）÷（計算値）−１）はΔ_b^’／Δ_spec＝１．３で約２％になる。
・ρ₂^’、ω_e2^*’はほぼ同値でΔ_b^’／Δ_specにほぼ比例する。

【0156】

図８の読取から、仕様値に対するパラメータ不確かさであるΔ_b’ ／Δ_spec、に対応するρ₂^’、ω_e2^*’の予測値を次式で近似することができる。

【0157】

【数81】

ここで近似誤差はΔ_b^’／Δ_spec＝１．３のときρ~₂^’が約０．８％、ω~_e2^*’が約１．８％になる。

【0158】

このように、Δ_bの増加によるζ_e2^*の劣化、ω_e2^*の変化およびρ₂の増加は、Δ_b’／Δ_specによる予測値から把握できることが分かる。

【0159】

２．３．２ ３次推定器
３次推定器は２次推定器に、舵角オフセット成分δ_oの推定を加えたもので、５次特性多項式は、（３４）式から、

【0160】

【数82】

になる。ここで、Ｂ₃（ｓ）は、（２８）式より、

【0161】

【数83】

を示す。

【0162】

２次推定器との対比を行って、２次推定器からの変化の特性を把握することによって、３次推定器の解析を行う。まず、（６５）式のパラメータ不確かさの多項式を、

【0163】

【数84】

の因数分解した形に置き換える。ここで、ζ_f ＝ ζ_e を用いて

【0164】

【数85】

を示す。２次推定器で用いた具体的数値とρ₂＝４．５８、ρ_o＝０．１から４ａ_zｃ_z／ｂ_z²＝０．３０５となり、４ａ_zｃ_z／ｂ_z²は小さいことが推定される。よって、

【0165】

【数86】

なる近似を用いると（√（１＋ε）≒１＋ε／２ （０＜｜ε｜＜１）））、Ｚ（ｓ）の根として次式を得る（ｚ₂＜ｚ₁とする）。

【0166】

【数87】

【0167】

【数88】

【0168】

（６８）式において分母は１に近いため、ｚ₁≒ω_eoであり、しかも、Ｄ_eo3^＊（ｓ）の特性根−ω_eo3^*は極−ω_eoからゼロ点ｚ₁の間にあるので、Ｄ_eo3^*（ｓ）≒Ｄ_eo（ｓ）と近似することができる。これより（６５）式は次式になる。

【0169】

【数89】

【0170】

【数90】

【0171】

（７１）式を見ると、Ｄ_c43^Δb（ｓ）は２次推定器の（５８）式のＤ_c4^Δb（ｓ）と同じ形になっていることが分かる。そこで、両者のΔ_bの係数を比較するために、（５９）式のｓ¹ 係数ｂ₁₂、ｓ⁰ 係数ｂ₀₂と（７１）式のｓ¹ 係数ａ_z、ｓ⁰ 係数ａ_zｚ₂とをそれぞれ比較すると、

【0172】

【数91】

【0173】

【数92】

【0174】

になる。（７２）式及び（７３）式の比がそれぞれ１より大きいことから、パラメータ不確かさΔ_bの係数は、３次推定器の５次特性多項式の方が、２次推定器の４次特性多項式に比べ、大きくなることが分かり、２次推定器におけるパラメータ不確かさΔ_bが等価的に大きくなったものとみなすことができる。そのため３次推定係数ρ₃は、ρ₂に対してパラメータ不確かさΔ_bの等価的増加分を修正する必要がある。Δ_bの係数比を具体的な数値を代入して計算すると（ρ_o＝０．１）

【0175】

【数93】

になる。上式から大きい値のｓ¹ 係数を選ぶこととすると、２次推定器の４次特性多項式においてパラメータ不確かさがａ_z /ｂ₁₂倍増加したことと等価になる。よって、前節の２次推定器の特性から（６４）式を用いて、ρ’₃の予測値を換算すると

【0176】

【数94】

と定めることができる。これよりρ₃、ω_e3^*の予測値を次式に定めることができる。

【0177】

【数95】

【0178】

２．３．３ ４次推定器
４次推定器は２次推定器に波浪成分ψ_wの推定を加えたもので、６次特性多項式は（３５）式から、

【0179】

【数96】

になる。ここで、（２９）式より

【0180】

【数97】

である。

【0181】

３次推定器での手法と同様に、２次推定器との対比を行って、４次推定器の解析を行う。そのため、Ｂ₄の係数ｂ₁₄と、Ｂ₂の係数ｂ₁₂との比較をすると、（７９）式から

【0182】

【数98】

となる。詳細は省略するが、ｆ₁Ｔ_snω_e2²／ｂ₁₂＞１であるため、上記式は、ｋ_wに関して右下がりの直線となる。

【0183】

波浪パラメータの中心角周波数比を

【0184】

【数99】

と定め、ζ_wn＝０．１として数値計算すると

を得る。最悪ケースのｂ₁₄ /ｂ₁₂＝１．１２１は（７４）式のａ_z /ｂ₁₂＝１．１３７とほぼ同等であるが、ｒ_w≫１でｂ₁₄はｂ₁₂に収束する。

【0185】

（７８）式でＤ_w(s)とＤ_ew(s)とは次の条件のとき

【0186】

【数100】

に近似できるので、２次推定器の特性と等価になり、４次推定係数ρ₄≒ρ₂，ω_e4^*≒ω_e2^*になる。上記の近似ができない場合、ρ₄は減衰係数差（ζ_e−ζ_wn）に比例し、周波数比ｒ_wに反比例する傾向を持つとし、その予測値を次式に定める。

【0187】

【数101】

【0188】

２．３．４ ５次推定器
５次推定器は２次推定器に舵角オフセット成分δ_o、波浪成分ψ_wの推定を加えたもので、３次及び４次の推定根特性を踏襲するものとし、７次特性多項式は（３６）式から、次式となる。

【0189】

【数102】

【0190】

上式の予測値を３次、４次の推定根特性から次式に定める。

【0191】

【数103】

【0192】

２．４ 波浪パラメータ条件
Ｄ_c7^Δb（ｓ）において、仕様を満たすがω_e5 ＝ ρ₅ω_f＞ω_wn になった場合、波浪パラメータ条件（表１）に反する。そのときζ_wn=ζ_e に変更して再度ρ₅ を求める。ノッチフィルター効果（外乱除去比ζ_wn /ζ_ew）を失うが、ρ₅ が小さくなるのでω_e5^* を下げローパスフィルタ効果を得ることができる。すなわち（８５）式で

【0193】

【数104】

になり、Ｄ_c5^Δb（ｓ）の３次推定器の特性に変更される。

【0194】

２．５ まとめ
２．４節の結果をまとめると、
・推定根のζ_e2^*, ω_e2^* はζ_spec付近で、ρ₂ にほぼ比例する。
・外乱モデルを含むρ₅とω_e5^* との予測値は、外乱モデルを含まないρ₂とω_e2^* によって求められる。

【0195】

７次特性根からρ₅ を計算するとき、ρ₅ とζ_e5^*, ω_e5^* との初期値が適切であると求解が確実になる。

【0196】

以下、その求解算法について説明する。

【0197】

３．求解算法
推定係数ρを求める算法を示す。求められたρはパラメータ不確かさを印加した閉ループ系において、仕様Δ_spec, ζ_specを満足する。算法は確実な求解法、簡単な処理と少ない計算量とのため、前章で導き出した性質を用いて、３つの段階により構築する。

【0198】

３．１ メインフロー
次の流れを採用する（図９）
（１）（５８）式でζ_e2^*＝ζ_spec になるρ₂，ω_e2^* を求める（Ｓ１０）。
（２）（８５）式でζ_e5^*＝ζ_spec になるρ₅，ω_e5^* を求めるために、（８６）式の予測値を初期値に設定する（Ｓ２０）。またω_e5＞ω_wn のとき（Ｓ３０でｙｅｓ）、ζ_wn＝ζ_e にして（Ｓ４０）、（２）を繰り返す。
（３）求めたρ₅をρ（＝ω_e／ω_f）として制御系を構成する（Ｓ５０）。

【0199】

３．２ ステップ２ サブフロー１
ρ（ρ₂、ρ₅）はその性質を利用して２段階で求める（図１０参照）。

【0200】

（１）ζ_spec を挟む範囲のρを求める：ζ_e^* はρに対して単調増加になる性質を利用する。初期値は

【0201】

【数105】

とおく。ステップ３に初期値を送るとζ_e^*，ω_e^* が戻される。反復ループ（添え字ｉ：回数）で

【0202】

【数106】

を更新する。上記にてｅ_ζi＝ζ_ei^*−ζ_spec，ｅ_ζi＞０のときブレークし、下限値と上限値とを設定する（図１０参照）。

【0203】

（２）次に、ζ_e^*＝ζ_spec となるρを求める。ζ_e^* はζ_spec 付近でρにほぼ直線近似になる性質を利用する。反復ループ（添え字ｉ：回数）で

【0204】

【数107】

を更新する。上記にて｜ｅ_ζi｜＜ε（ε：微小項）のときブレークする。｜ｅ_ζi｜≧εのときに、範囲を次式に置き換えて反複ループを実行する（図１０参照）。

【0205】

【数108】

【0206】

３．３ ステップ３
初期値を用いて次の処理を行い、ζ_e^*とω_e^*とを取り出す。
（１）４次式は代数解法によって、最小の減衰係数に対応する２次式の係数を求める。
（２）７次式は任意の数値解法によって求める。具体的には、ベアストウ法 と組立除法とによって、２次式の係数を求め、残りの５次式を因数分解する場合は ニュートン−ラプソン法（初期値ω_eo）と代数解法とを用いる。

【0207】

なお、上記の解法は公知のため説明を省く。

【0208】

４．計算結果
本発明を数値計算によって検証する。船体パラメータおよび制御パラメータを表３にまとめる。ただし設計パラメータについては表１を参照されたい。また、仕様は（４０）、（４２）式に従っている。

【0209】

【表3】

【0210】

４．１ ２次推定器
表３から以下のことが分かる。
・ρ₂の変化は微分ゲインｆ₂ のそれより小さいが、ω_fＴ_sn に関係する。ρ₂ は仕様から判断すると実用範囲の値である。

【0211】

４．２ ３次推定器
図１１を参照すると
・ａ_z /ｂ₁₂：ω_fＴ_sn ＝１ の場合は他の場合より１５〜３０％大きくなるが、ω_fＴ_sn ＞１ の場合はほぼ近い値を示す。
・（−ａ_zｚ₂）/ｂ₀₂ ：近似値はａ_z /ｂ₁₂ の９０％弱になり、ρ_o が大きくなると近似誤差が増加する。
・誤差ρ₃~’，誤差ω_e3~’：ρ_o≦０．２とすれば、両者とも５％以内の誤差に収まる。

【0212】

よって、３次推定器を２次推定器に近似し予測値を求める方法は妥当であり、予測値ρ₃~’，ω_e3~^*’ は有効であると言える。

【0213】

４．３ ４次推定器
図１２（ω_fＴ_sn ＝１．５８を利用）を参照すると
・ｂ₁₄ /ｂ₁₂ 、ｋ_w：共にｒ_w ＝１で１に近く、ｒ_w ＝１０でほぼ１に収斂する。
・ζ_e4^*’，ω_e4^*’：Ｂ₄（ｓ）のみの場合（波浪成分Ｄ_ew（ｓ），Ｄ_w（ｓ）を省く）は波浪成分を含む場合より変化量が小さい。

【0214】

図１３より、ρ₄^’ の予測値はζ_wn ＝０．１かつｒ_w が低減で誤差を生じるが、それ以外では誤差は小さいと見なせる。ω_e4^’ の予測値は一定のため、低い領域で誤差を生じる。

【0215】

これに対して、波浪成分による影響が支配的になり、予測値ρ₄^’ ，ω_e4^’ は実用精度をもつ。

【0216】

４．４ ５次推定器
・図１４は図１３にρ_o を追加した場合で、ほぼ同様の傾向を示し、オフセットは予測値とほぼ等価である。
・図１５はζ_wn＝０．１，ω_fＴ_sn を変化させた場合で、ω_fＴ_sn ＝１のｒ_w の低減で急増する。

【0217】

予測値ρ₅^’ ，ω_e5^’ はｒ_w の低減で波浪成分誤差をもつが、波浪パラメータ条件（破線より右側部分が条件を満足する領域）から予測値として実用精度をもつことが分かる。

【0218】

４．５ 計算時間の比較
本発明の求解算法の計算時間を標準算法と比較する。標準算法は既知の推定係数ρ_sol(
ｒ_w) から探索範囲［ρ_low，ρ_high］＝ρ_sol (ｒ_w)［1/ｃ_p，ｃ_p］を与えて（Matlab のfminbnd を利用）、（１２）式の特性行列から特性根を求め、最小の減衰係数を取り出し、それが仕様と一致する推定係数を求める。計算結果を表４にまとめる。条件：図１４，ζ_wn＝０．１である。

【0219】

【表4】

【0220】

本発明の計算時間は０．５８２，０．５８１，０．５８５[ｓ]，平均０．５８３[ｓ] になり、標準算法の計算時間に比べ、１桁程度速いことが確認できた。

【0221】

５．まとめ
推定器を含んだ保針ループに対して、推定根に影響を与えるパラメータ不確かさを特定し、仕様を定め、推定根特性を仕様付近で把握することで、ρの予測値を求めた。外乱モデルなしの２次推定根特性を基に外乱モデルありの推定根特性を解析した。求解算法は数値解と推定根特性とを利用し予測値を初期値に採用した。

【0222】

数値計算によって、推定根特性と予測値との妥当性を確認し、求解算法の有効性を標準算法の計算時間と比べ確認した。

【符号の説明】

【0223】

１０ 自動操舵装置
１２ フィードバック制御器
１４ 操舵機
１６ 船体
１８ 制御対象
２２ 状態フィードバック制御器
２４ 推定器

【図1】

【図1A】

【図2】

【図3A】

【図3B】

【図3C】

【図3D】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】