(57)【要約】【課題】従来に比して推力発生装置の発生推力を高精度に制御することができる推力配分方法及び推力配分装置を提供する。【解決手段】推力発生装置１〜３の非線形特性及び動特性を考慮した航走体１の推力に関する力及びモーメントのバランス関係を予め設定しておき、また、各推力発生装置の発生推力を可及的に小さくするための評価関数を規定する。前記バランス関係を拘束条件とし、前記評価関数の値を最小とするような各推力発生装置１〜３のプロペラ推力及び舵２２，３２の舵角並びに舵角速度を演算する。【選択図】　図１
2004042885.tif 000002
【特許請求の範囲】
【請求項１】
操作量に対応する推力を発生する複数の推力発生装置を備える航走体を、与えられた目標推力で航走又は定点保持させるために、夫々の推力発生装置に与える操作量を演算して、航走体の推力配分を行う推力配分方法において、
夫々の推力発生装置の非線形特性に基づいて予め定められた、夫々の推力発生装置が発生する推力及び航走体の推力に関する力及びモーメントのバランス関係と、
推力発生装置が発生する推力を可及的に小さくするための、夫々の推力発生装置の状態に関する状態量及び／又は発生推力に関する所定の評価関数とに基づいて、
前記目標推力との間で前記バランス関係を満たし、且つ前記評価関数に対して所定の条件を満たす夫々の推力発生装置の操作量を演算することを特徴とする推力配分方法。
【請求項２】
前記バランス関係は、夫々の推力発生装置の非線形特性及び動特性に基づいて予め定められていることを特徴とする請求項１に記載の推力配分方法。
【請求項３】
前記バランス関係は、（１）〜（３）式で与えられることを特徴とする請求項１又は２に記載の推力配分方法。
【数１】
2004042885.tif 000094
【数２】
2004042885.tif 000095
【数３】
2004042885.tif 000096
但し、Ｘ：航走体の推力の前後方向成分
Ｙ：航走体の推力の左右方向成分
Ｎ：航走体の旋回方向のモーメント
ｆｘｉ：ｉ番目の推力発生装置の発生推力の前後方向成分
ｆｙｉ：ｉ番目の推力発生装置の発生推力の左右方向成分
ｋ１ｉ，ｋ２ｉ：係数
ｉ：１〜ｎの自然数
【請求項４】
夫々の推力発生装置の発生推力Ｆ（ｔ）＝［ｆｘ１　…　ｆｘｎｆｙ１　…　ｆｙｎ］Ｔは、夫々の推力発生装置の操作量をｕ（ｔ）＝［ｕ１　…　ｕｍ］Ｔ（ｍは正の整数）としたときに、（４）式の非線形方程式で与えられることを特徴とする請求項３に記載の推力配分方法。
【数４】
2004042885.tif 000097
但し、ｔ：時間変数
ｘ：（５）式の非線形微分方程式で与えられる状態量
【数５】
2004042885.tif 000098
【請求項５】
前記操作量ｕは、（６）式で与えられる関係を満たすことを特徴とする請求項４に記載の推力配分方法。
【数６】
2004042885.tif 000099
但し、ｕｊｍａｘ：定数
ｊ：１〜ｍの自然数
【請求項６】
前記状態量ｘは、（７）式で与えられる関係を満たすことを特徴とする請求項４又は５に記載の推力配分方法。
【数７】
2004042885.tif 000100
但し、ｘｉｍａｘ：定数
【請求項７】
前記推力発生装置が旋回式スラスタであり、前記バランス関係は、（８）〜（１０）式で与えられることを特徴とする請求項４に記載の推力配分方法。
【数８】
2004042885.tif 000101
【数９】
2004042885.tif 000102
【数１０】
2004042885.tif 000103
但し、Ｔｉ：ｉ番目の推力発生装置のプロペラ推力
θｉ：ｉ番目の推力発生装置の旋回角度
ｘｏｉ：ｉ番目の推力発生装置の航走体重心に対する前後方向の位置
ｙｏｉ：ｉ番目の推力発生装置の航走体重心に対する左右方向の位置
【請求項８】
前記評価関数Ｊは、（１１）式で与えられ、前記所定の条件は、前記評価関数Ｊの値を最小化又は最大化することとして与えられることを特徴とする請求項４乃至７の何れかに記載の推力配分方法。
【数１１】
2004042885.tif 000104
但し、φ，Ｌ：所定階以上微分可能なスカラー値関数
Ｔ：定数
ｔ’：ｔ≦ｔ’≦ｔ＋Ｔを満たす変数
【請求項９】
前記評価関数Ｊは、状態量ｘ及び／又は操作量ｕに関する重み係数を含み、
該重み係数をオンラインで変更することを特徴とする請求項８に記載の推力配分方法。
【請求項１０】
推力発生装置に故障が発生した場合に、故障が発生した推力発生装置の状態量ｘ及び／又は操作量ｕに関する重み係数を変更することを特徴とする請求項９に記載の推力配分方法。
【請求項１１】
操作量に対応する推力を発生する複数の推力発生装置を備える航走体を、与えられた目標推力で航走又は定点保持させるために、夫々の推力発生装置に与える操作量を演算して、航走体の推力配分を行う推力配分装置において、
夫々の推力発生装置の非線形特性に基づいて予め定められた、夫々の推力発生装置が発生する推力及び航走体の推力に関する力及びモーメントのバランス関係と、
推力発生装置が発生する推力を可及的に小さくするための、夫々の推力発生装置の状態に関する状態量及び／又は発生推力に関する所定の評価関数とに基づいて、
前記目標推力との間で前記バランス関係を満たし、且つ前記評価関数に対して所定の条件を満たす夫々の推力発生装置の操作量を演算する演算手段を備えることを特徴とする推力配分装置。
【請求項１２】
前記バランス関係は、夫々の推力発生装置の非線形特性及び動特性に基づいて予め定められていることを特徴とする請求項１１に記載の推力配分装置。
【請求項１３】
前記バランス関係は、（１）〜（３）式で与えられることを特徴とする請求項１１又は１２に記載の推力配分装置。
【数１２】
2004042885.tif 000105
【数１３】
2004042885.tif 000106
【数１４】
2004042885.tif 000107
但し、Ｘ：航走体の推力の前後方向成分
Ｙ：航走体の推力の左右方向成分
Ｎ：航走体の推力の旋回方向成分
ｆｘｉ：ｉ番目の推力発生装置の発生推力の前後方向成分
ｆｙｉ：ｉ番目の推力発生装置の発生推力の左右方向成分
ｋ１ｉ，ｋ２ｉ：係数
ｉ：１〜ｎの自然数
【請求項１４】
夫々の推力発生装置の発生推力Ｆ（ｔ）＝［ｆｘ１　…　ｆｘｎ　ｆｙ１　…　ｆｙｎ］Ｔは、夫々の推力発生装置の操作量をｕ（ｔ）＝［ｕ１　…ｕｍ］Ｔ（ｍは正の整数）としたときに、（４）式の非線形方程式で与えられることを特徴とする請求項１３に記載の推力配分装置。
【数１５】
2004042885.tif 000108
但し、ｔ：時間変数
ｘ：（５）式の非線形微分方程式で与えられる状態量
【数１６】
2004042885.tif 000109
【請求項１５】
前記操作量ｕは、（６）式で与えられる関係を満たすことを特徴とする請求項１４に記載の推力配分装置。
【数１７】
2004042885.tif 000110
但し、ｕｊｍａｘ：定数
ｊ：１〜ｍの自然数
【請求項１６】
前記状態量ｘは、（７）式で与えられる関係を満たすことを特徴とする請求項１４又は１５に記載の推力配分装置。
【数１８】
2004042885.tif 000111
但し、ｘｉｍａｘ：定数
【請求項１７】
前記推力発生装置が旋回式スラスタであり、前記バランス関係は、（８）〜（１０）式で与えられることを特徴とする請求項１４に記載の推力配分方法。
【数１９】
2004042885.tif 000112
【数２０】
2004042885.tif 000113
【数２１】
2004042885.tif 000114
但し、Ｔｉ：ｉ番目の推力発生装置のプロペラ推力
θｉ：ｉ番目の推力発生装置の旋回角度
ｘｏｉ：ｉ番目の推力発生装置の航走体重心に対する前後方向の位置
ｙｏｉ：ｉ番目の推力発生装置の航走体重心に対する左右方向の位置
【請求項１８】
前記評価関数Ｊは、（１１）式で与えられ、前記所定の条件は、前記評価関数Ｊの値を最小化又は最大化することとして与えられることを特徴とする請求項１４乃至１７の何れかに記載の推力配分装置。
【数２２】
2004042885.tif 000115
但し、φ，Ｌ：所定階以上微分可能なスカラー値関数
Ｔ：定数
ｔ’：ｔ≦ｔ’≦ｔ＋Ｔを満たす変数
【請求項１９】
前記評価関数Ｊは、状態量ｘ及び／又は操作量ｕに関する重み係数を含み、
該重み係数をオンラインで変更可能な変更手段を更に備えることを特徴とする請求項１８に記載の推力配分装置。
【請求項２０】
推力発生装置の故障の発生を検出する検出手段を更に備え、
前記変更手段は、該検出手段が推力発生装置の故障の発生を検出した場合に、故障が発生した推力発生装置の状態量ｘ及び／又は操作量ｕに関する重み係数を変更すべくなしてあることを特徴とする請求項１９に記載の推力配分装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、船舶及び飛行船等の航走体が所要の推力を発生するために、該航走体に設けられた複数の推力発生装置に与えるべき操作量を演算して、航走体の推力配分を行う推力配分方法及び該推力配分方法の実施に直接使用する推力配分装置に関する。
【０００２】
【従来の技術】
カーフェリー及び潜水艇等の船舶並びに飛行船等の航走体には、複数の推力発生装置（スラスタ）が設けられていることが多い。かかる航走体を航走又は定点保持させるためには、複数の推力発生装置を駆動する必要があるが、航走体全体で一の推力を発生させるのに、夫々の推力発生装置の発生推力の組み合わせが複数組存在する場合がある。従って、推力発生装置の消費燃料を可及的に少なくするために、如何に各推力発生装置の発生推力の配分を行うかが重要となっている。このような要請から、航走体全体の目標とする推力を発生させるために、夫々の推力発生装置が発生する推力を配分する種々の方法が提案されている。
【０００３】
この種の方法の一つとして、特許文献１に開示されたものがある。この推力配分方法は、夫々の推力発生装置の発生推力と、航走体全体の推力とに関する力及びモーメントのバランス関係（拘束条件）を満たしつつ、夫々の推力発生装置の発生推力に関する所定の評価関数の値を最小化するような夫々の推力発生装置の発生推力を演算するものである。
【０００４】
また、特許文献２には、航走体が首振り式の推力発生装置を備える場合に、夫々の推力発生装置の定格推力及び推力発生装置の動作角度範囲（首振り角範囲）の制限を満たしつつ、消費エネルギを最小にするような夫々の推力発生装置の発生推力を演算する推力配分方法が開示されている。この推力配分方法は、力及びモーメントのバランス関係を満たしつつ、各推力発生装置の発生推力に関する評価関数の値を最小化するような夫々の推力発生装置の発生推力を演算し、演算結果の推力が定格推力を越えている場合には、当該推力発生装置に対応するパラメータを変更して再度推力を演算し、また推力発生装置にその方位角が禁止角度範囲内にあるものが存在する場合には、当該推力発生装置の方位角を禁止角度範囲外に設定するものである。
【０００５】
また、特許文献３には、航走体が首振り式の推力発生装置を備える場合に、夫々の推力発生装置の首振り角度が最小となるように、夫々の推力発生装置の発生推力を演算する推力配分方法が開示されている。この推力配分方法は、力及びモーメントのバランス関係を満たしつつ、各推力発生装置の首振り角度に関する評価関数の値を最小化することにより、夫々の推力発生装置の首振り角度を可及的に小さくするような発生推力を演算するものである。
【０００６】
【特許文献１】
特開平１−３１１９９７号公報
【特許文献２】
特開２００１−２１９８９９号公報
【特許文献３】
特開平９−２３９３号公報
【０００７】
【発明が解決しようとする課題】
ところが、上述の如き従来の推力配分方法にあっては、拘束条件たる力及びモーメントのバランス関係が、推力発生装置の非線形特性や動特性に基づいて定められていなかった。一般的に、推力発生装置は非線形特性を示す場合が多い。これは例えば、推力発生装置がプロペラである場合に、プロペラの回転速度と発生推力との関係が非線形であるようなものが考えられる。このような場合には、演算によって得られた推力を発生するために、単純にこの推力に比例した回転速度（操作量）でプロペラを駆動したときには、目標の推力に対して実際に発生した推力に大きな誤差が生じることとなる。このため、従来の推力配分方法では、推力発生装置の発生推力を高精度に制御することができなかった。
【０００８】
また、推力発生装置の特性が遅れ要素等の動特性を含む場合もある。これは例えば、比較的応答速度が低い舵を備える推力発生装置等が考えられる。このような場合には、演算によって得られた推力を発生するために必要な舵角に変角するように舵を駆動したときには、所要の舵角まで舵が変角するのに時間がかかることとなる。このため、このような場合にも、従来の推力配分方法では、推力発生装置の発生推力を高精度に制御することができなかった。
【０００９】
また、上述した評価関数には、推力発生装置の状態を示す状態量及び／又は操作量に掛け合わせる重み係数を含むことが多い。推力発生装置の一部が故障した場合には、この推力発生装置に関する重み係数を変更すること等によって、故障が発生した推力発生装置を動作させないように設定し、残りの推力発生装置のみで推力配分を行うようにしなければ、航走体全体での目標推力を正確に得ることができない。ところが、従来の推力配分方法には、この重み係数をオンラインで変更するものがなかったため、推力発生装置が故障した場合に、残りの推力発生装置の発生推力を正確に求めることができず、実質的に推力配分が不可能であった。
【００１０】
本発明は、斯かる事情に鑑みてなされたものであり、従来に比して推力発生装置の発生推力を高精度に制御することができる推力配分方法及び該推力配分方法の実施に直接使用する推力配分装置を提供することを目的とする。
【００１１】
また、本発明の他の目的は、推力発生装置の一部が故障した場合であっても、残りの推力発生装置の発生推力のみを用いた正確な演算をすることができる推力配分方法及び該推力配分方法に直接使用する推力配分装置を提供することにある。
【００１２】
【課題を解決するための手段】
上記課題を解決するため、本発明に係る推力配分方法は、操作量に対応する推力を発生する複数の推力発生装置を備える航走体を、与えられた目標推力で航走又は定点保持させるために、夫々の推力発生装置に与える操作量を演算して、航走体の推力配分を行う推力配分方法において、夫々の推力発生装置の非線形特性に基づいて予め定められた、夫々の推力発生装置が発生する推力及び航走体の推力に関する力及びモーメントのバランス関係と、推力発生装置が発生する推力を可及的に小さくするための、夫々の推力発生装置の状態に関する状態量及び／又は発生推力に関する所定の評価関数とに基づいて、前記目標推力との間で前記バランス関係を満たし、且つ前記評価関数に対して所定の条件を満たす夫々の推力発生装置の操作量を演算することを特徴とする。
【００１３】
また、本発明に係る推力配分装置は、操作量に対応する推力を発生する複数の推力発生装置を備える航走体を、与えられた目標推力で航走又は定点保持させるために、夫々の推力発生装置に与える操作量を演算して、航走体の推力配分を行う推力配分装置において、夫々の推力発生装置の非線形特性に基づいて予め定められた、夫々の推力発生装置が発生する推力及び航走体の推力に関する力及びモーメントのバランス関係と、推力発生装置が発生する推力を可及的に小さくするための、夫々の推力発生装置の状態に関する状態量及び／又は発生推力に関する所定の評価関数とに基づいて、前記目標推力との間で前記バランス関係を満たし、且つ前記評価関数に対して所定の条件を満たす夫々の推力発生装置の操作量を演算する演算手段を備えることを特徴とする。
【００１４】
上記発明によれば、夫々の推力発生装置が発生する推力及び航走体の推力に関する力及びモーメントのバランス関係を、夫々の推力発生装置の非線形特性に基づいて予め定めておき、航走体の目標推力に対して、この力及びモーメントのバランス関係を満たすような、夫々の推力発生装置の発生推力の中から、評価関数に対して所定の条件を満たす推力を求めることとなる。これにより、操作量と発生推力との関係が線形性を示さない推力発生装置を航走体が備えている場合であっても、目的とする発生推力を得るために必要な操作量を正確に演算することができ、推力発生装置の発生推力を従来に比して高精度に制御することが可能となる。
【００１５】
また、上記発明においては、前記バランス関係を、夫々の推力発生装置の非線形特性及び動特性に基づいて予め定められたものとすることが望ましい。
【００１６】
この場合には、操作量と発生推力との関係が線形性を示さず、しかも動作指示が与えられてから指示された推力を発生するまでに時間的な遅れを有する、即ち動特性を有する推力発生装置を航走体が備えている場合であっても、目的とする発生推力を得るために必要な操作量を正確に演算することができ、推力発生装置の発生推力を従来に比して高精度に制御することが可能となる。
【００１７】
また、上記発明においては、前記バランス関係を、（１）〜（３）式で与えることもできる。
【００１８】
【数２３】
2004042885.tif 000003
【数２４】
2004042885.tif 000004
【数２５】
2004042885.tif 000005
但し、Ｘ：航走体の推力の前後方向成分
Ｙ：航走体の推力の左右方向成分
Ｎ：航走体の旋回方向のモーメント
ｆｘｉ：ｉ番目の推力発生装置の発生推力の前後方向成分
ｆｙｉ：ｉ番目の推力発生装置の発生推力の左右方向成分
ｋ１ｉ，ｋ２ｉ：係数
ｉ：１〜ｎの自然数
【００１９】
また、この場合において、夫々の推力発生装置の発生推力Ｆ（ｔ）＝［ｆｘ１　…　ｆｘｎ　ｆｙ１　…　ｆｙｎ］Ｔを、夫々の推力発生装置の操作量をｕ（ｔ）＝［ｕ１　…　ｕｍ］Ｔ（ｍは正の整数）としたときに、（４）式の非線形方程式で与えることもできる。
【００２０】
【数２６】
2004042885.tif 000006
但し、ｔ：時間変数
ｘ：（５）式の非線形微分方程式で与えられる状態量
【数２７】
2004042885.tif 000007
【００２１】
これにより、非線形性を示す推力発生装置の動特性を従来に比して正確にモデル化することができ、夫々の推力発生装置の発生推力をより高精度に制御することが可能となる。
【００２２】
また、この場合においては、前記操作量ｕを、（６）式で与えられる関係を満たすように定めることもできる。
【００２３】
【数２８】
2004042885.tif 000008
但し、ｕｊｍａｘ：定数
ｊ：１〜ｍの自然数
【００２４】
また、上記発明においては、前記状態量ｘを、（７）式で与えられる関係を満たすように定めることもできる。
【００２５】
【数２９】
2004042885.tif 000009
但し、ｘｉｍａｘ：定数
【００２６】
例えば、１つのプロペラと該プロペラの後方に配置された１つの舵とによって構成された推力発生装置では、プロペラによって後方へ流体を送出する場合であって、プロペラによる流体の送出方向に交差する方向へ操舵されているときには、舵に揚力が発生し、これによって横方向の力（推力の横方向成分）が発生するが、プロペラが逆回転することによって流体を前方へ送出する場合であって、流体の送出方向に交差する方向に操舵されているときには、送出される流体の流れによる舵揚力は殆ど発生せず、よって横方向への力は殆ど発生しない。このような非線形モデルは、前述した（４）式によって表現することができる。従って、推力発生装置の非線形特性を従来に比して正確にモデル化することができ、夫々の推力発生装置の発生推力をより高精度に制御することが可能となる。
【００２７】
また、上記発明においては、前記推力発生装置を旋回式スラスタとし、前記バランス関係を、（８）〜（１０）式で与えることもできる。
【００２８】
【数３０】
2004042885.tif 000010
【数３１】
2004042885.tif 000011
【数３２】
2004042885.tif 000012
但し、Ｔｉ：ｉ番目の推力発生装置のプロペラ推力
θｉ：ｉ番目の推力発生装置の旋回角度
ｘｏｉ：ｉ番目の推力発生装置の航走体重心に対する前後方向の位置
ｙｏｉ：ｉ番目の推力発生装置の航走体重心に対する左右方向の位置
【００２９】
また、上記発明においては、前記評価関数Ｊを、（１１）式で与え、前記所定の条件を、前記評価関数Ｊの値を最小化又は最大化することとして与えることもできる。
【００３０】
【数３３】
2004042885.tif 000013
但し、φ，Ｌ：所定階以上微分可能なスカラー値関数
Ｔ：定数
ｔ’：ｔ≦ｔ’≦ｔ＋Ｔを満たす変数
【００３１】
このように、評価区間が移動する評価関数を最小化するフィードバック制御は、Ｒｅｃｅｄｉｎｇ　Ｈｏｒｉｚｏｎ制御と呼ばれている。
ここで、現在時刻ｔに対して常にτ＝０が対応する仮想的な時間変数τを導入する。これにより、上述の如き問題は、（１２）〜（１４）式で表すような時刻ｔに関するτ軸上の有限時間最適制御問題に帰着される。
【００３２】
【数３４】
2004042885.tif 000014
【数３５】
2004042885.tif 000015
【数３６】
2004042885.tif 000016
【００３３】
但し、ここでは下付き添字は偏微分を表すものとしており、また仮想的時間軸τ上の変数に添字＊を付加するものとしている。そして、ｘ＊（τ，ｔ），ｕ＊（τ，ｔ）が、ｘ（ｔ）から出発してＴだけ未来までの最適軌道及び最適制御入力である。また、実際の制御入力は、（１５）式で与えられる。
【００３４】
【数３７】
2004042885.tif 000017
【００３５】
仮想的時間軸τに沿っては、通常の最適制御問題になる。ここで、λ＊，ρ＊を未定乗数として、ハミルトニアンＨを（１６）式で定義する。
【００３６】
【数３８】
2004042885.tif 000018
但し、Ｃ：（１７）式で与えられるベクトル値関数
【数３９】
2004042885.tif 000019
【００３７】
これによって、（１８）〜（２３）式のようなオイラー・ラグランジュ方程式が得られる。
【００３８】
【数４０】
2004042885.tif 000020
【数４１】
2004042885.tif 000021
【数４２】
2004042885.tif 000022
【数４３】
2004042885.tif 000023
【数４４】
2004042885.tif 000024
【数４５】
2004042885.tif 000025
【００３９】
以上の如きＲｅｃｅｄｉｎｇ　Ｈｏｒｉｚｏｎ制御については、例えば、計測と制御第３６巻第１１号（１９９７年１１月号）の７７６〜７８３頁等に、実時間で解くことが可能な解法が提案されているので、このような解法を用いて、本発明に係る推力配分方法及び推力配分装置を構成すれば、オンライン、即ち実時間内で航走体の推力配分を行うことが可能となる。
【００４０】
また、上記発明においては、前記評価関数Ｊを、状態量ｘ及び／又は操作量ｕに関する重み係数を含むように規定しておき、該重み係数をオンラインで変更するように構成することが望ましい。
【００４１】
これにより、航走体が備える推力発生装置の一部に故障が発生した場合に、故障が発生した推力発生装置に係る状態量ｘ及び／又は操作量ｕに関する重み係数を適宜変更することで、故障が発生した推力発生装置を動作させないように設定し、残りの推力発生装置で推力配分を継続して行うことができる。
【００４２】
また、この場合には、推力発生装置の故障の発生を検出するように構成し、故障の発生が検出された推力発生装置の状態量ｘ及び／又は操作量ｕに関する重み係数を変更するように構成することが望ましい。これにより、推力発生装置に故障が発生した場合の重み係数の変更を自動的に行わせることができる。
【００４３】
【発明の実施の形態】
以下、本発明の実施の形態に係る推力配分装置について、図面を参照しながら具体的に説明する。
【００４４】
（実施の形態１）
図１は、本発明の実施の形態１に係る推力配分装置が搭載された船舶の構成を示す模式図である。図１に示すように、船舶１０は、３つの推力発生装置１，２，３を備えている。このうち、推力発生装置１は、船舶１０の船首側の部分に設けられた所謂バウスラスタであり、船体の横方向（図中Ｙ方向）へ推力を発生することができるようになっている。
【００４５】
また、推力発生装置２，３は、夫々プロペラ２１，３１と、これらのプロペラ２１，３１の後方に設けられた舵２２，３２とを有している。推力発生装置２は、船尾の中央から左側の部分に設けられており、推力発生装置３は、船尾の中央から右側の部分に設けられている。プロペラ２１，３１は、夫々船体の前後方向へ延びた回転軸を中心に回転することが可能であり、これらが回転することによって、船体の前後方向（図中Ｘ方向）へ推力を発生することが可能である。また、舵２２，３２は、船体の左右方向へ揺動することが可能である。プロペラ２１，３１によって後方へ送られた海水等の流体の流れ方向に対し、舵２２，３２を傾斜させることにより、その主面（舵面）に揚力が発生し、これによって推力発生装置２，３からは、舵２２，３２の傾斜角（舵角）の直交方向へ推力が発生することとなる。
【００４６】
図２は、本発明の実施の形態１に係る推力配分装置が搭載された船舶１０の構成を示すブロック図である。図２に示すように、船舶１０は、推力発生装置１〜３の他に、駆動装置４１，４２，推力配分装置５，検出装置６，設定装置７，目標推力演算装置８及び制御装置９を備えている。
【００４７】
駆動装置４１は、推力発生装置１に接続されている。この駆動装置４１は、例えばエンジン及び変速機等から構成されており、推力発生装置１を所要の回転数で回転させるようになっている。
【００４８】
また、駆動装置４２は、推力発生装置２，３に夫々接続されている。この駆動装置４２は、例えばエンジン、変速機、差動装置、及び油圧装置等から構成されており、プロペラ２１，３１を、夫々所要の回転数によって回転させ、舵２２，３２を、夫々所要の舵角まで傾斜させるようになっている。また、駆動装置４２は、プロペラ２１，３１を相異なる回転数で回転させることが可能であり、舵２２，３２を相異なる舵角まで傾斜させることが可能である。
【００４９】
検出装置６は、船舶１０の現在位置を検出するＧＰＳ等の位置検出手段６１と、船舶１０の船速（対地速度）を検出する速度検出手段６２と、船舶１０の船首の方位角を検出するジャイロコンパス等の方位角検出手段６３とを備えている。また、検出装置６は、目標推力演算装置８に接続されており、検出した船舶１０の現在位置、船速及び船首の方位角の各情報を目標推力演算装置８に出力するようになっている。
【００５０】
設定装置７は、船舶１０の位置、船速及び船首の方位角を設定するようになっている。かかる設定装置７は、図示しない操作装置を具備しており、この操作装置を船員等が操作することによって、船員等が船舶１０の位置、船速及び船首の方位角を設定するようになっている。なお、船員等が操作装置を操作することにより、船舶１０の位置、船速及び船首の方位角を設定する構成に限定されるものではなく、例えば、設定装置７が自動的に船舶１０の位置、船速及び船首の方位角を設定する構成としてもよい。
【００５１】
このような設定装置７は、目標推力演算装置８に接続されており、設定した船舶１０の位置、船速及び船首の方位角の各情報を目標推力演算装置８に出力するようになっている。
【００５２】
目標推力演算装置８は、ＣＰＵ、ＲＯＭ及びＲＡＭ等から構成されており、検出装置６から入力した船舶１０の現在位置、船速及び船首の方位角の各情報並びに設定装置７から入力した船舶１０の位置、船速及び船首の方位角の各情報に基づいて、船舶１０全体の目標推力の前後方向成分及び左右方向成分並びに船舶１０の旋回方向の目標モーメントを演算する。このような目標推力演算装置８は、推力配分装置５に接続されており、演算結果の目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントを推力配分装置５に出力するようになっている。
【００５３】
なお、このような目標推力演算装置８による演算結果を推力配分装置５に入力する構成に限定されるものではなく、例えば、船員等がジョイスティック等の入力装置を用いて目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントを直接推力配分装置５に入力する構成であってもよい。
【００５４】
推力配分装置５は、後述するように、目標推力演算装置８から入力した目標推力の前後方法成分及び左右方向成分並びに旋回方向の目標モーメントに基づいて、推力発生装置１〜３の夫々のプロペラ推力及び舵２２，３２の夫々の舵角を演算する。このような推力配分装置５は、制御装置９に接続されており、演算結果の各プロペラ推力及び各舵角を制御装置９に出力するようになっている。
【００５５】
制御装置９は、駆動装置４１，４２に接続されている。また、船舶１０には、推力発生装置１のプロペラ及びプロペラ２１，３１の回転数を夫々検出する回転数センサと、舵２２，３２の舵角を夫々検出する舵角センサとが設けられている。これらの回転数センサ及び舵角センサは、制御装置９に接続されており、検出結果の回転数及び舵角の情報を、夫々制御装置９へ出力するようになっている。
【００５６】
また、制御装置９は、ＣＰＵ、ＲＯＭ及びＲＡＭ等から構成されており、推力配分装置５から入力した各プロペラ推力及び各舵角と、各種センサから入力した現在の各プロペラの回転数及び各舵角とに基づいて、推力発生装置１〜３の各プロペラ推力及び各舵角が、推力配分装置５から出力された各プロペラ推力及び各舵角に追従するように、例えばＰＩＤ制御等を用いて、推力発生装置１〜３の各プロペラ推力及び各舵角を制御する。この制御装置９から駆動装置４１，４２に対して制御信号が出力され、推力発生装置１〜３が夫々駆動される。
【００５７】
更に、制御装置９は、各種センサの測定値から得られた推力発生装置１〜３のプロペラ推力及び舵角の内の何れかが、推力配分装置５から入力したプロペラ推力及び舵角に追従しているか否かを判別し、追従していない場合には、その推力発生装置が故障したものとみなして、故障発生の情報を推力配分装置５へ出力するようにもなっている。
【００５８】
次に、推力配分装置５の構成を更に詳しく説明する。推力配分装置５は、ＣＰＵ５１，ＲＯＭ５２，ＲＡＭ５３，及び入出力インタフェース５４から構成されている。
【００５９】
ＲＯＭ５２には、コンピュータプログラム及びこのコンピュータプログラムで使用するデータ等が格納されている。ＣＰＵ５１は、ＲＯＭ５２に格納されたコンピュータプログラム又はＲＡＭ５３にロードされたコンピュータプログラムを実行することができる。また、ＲＡＭ５３は、ＣＰＵ５１がコンピュータプログラムを実行するときの作業領域として使用される。
【００６０】
ＣＰＵ５１，ＲＯＭ５２，ＲＡＭ５３，及び入出力インタフェース５４は、バスによって接続されており、相互にデータの転送が可能である。また、入出力インタフェース５４は、目標推力演算装置８及び制御装置９に接続されており、目標推力演算装置８及び制御装置９との間で各種のデータを送受信することが可能である。
【００６１】
図３は、本発明の実施の形態１に係る推力配分装置の演算処理の流れを示すフローチャートである。推力配分装置５は、ＣＰＵ５１がコンピュータプログラムを実行することにより、以下に説明するような処理を行うことができるようになっている。
【００６２】
まず、ＣＰＵ５１は、船舶１０の目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントを取得する（ステップＳ１）。この処理は、入出力インタフェース５４が受信した目標推力演算装置８の出力信号値を、ＣＰＵ５１が参照することによって行われる。
【００６３】
次に、ＣＰＵ５１は、推力発生装置１〜３の故障の発生の検出処理を行う（ステップＳ２）。この処理は、入出力インタフェース５４が制御装置９から受信した故障発生の情報を、ＣＰＵ５１が取得することによって行われる。ＣＰＵ５１は、この情報を取得した場合には、この情報によって特定される推力発生装置に故障が発生したものとみなす。
【００６４】
ステップＳ２において、推力発生装置１〜３の何れにも故障の発生が検出されなかった場合、ＣＰＵ５１は、ステップＳ４へ処理を移す。一方、ステップＳ２において、推力発生装置１〜３の何れかの故障の発生が検出された場合、ＣＰＵ５１は、重み係数の変更処理を行う（ステップＳ３）。この処理は、ＣＰＵ５１が、後述する（３０）式における重み係数ｓｆ１，ｓｆ２，ｒ１〜ｒ５，ｑ１，ｑ２のうち、故障が発生した推力発生装置に係るものの値を、大幅に増加させることによってなされる。
【００６５】
また、ＣＰＵ５１は、ステップＳ１にて取得した目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントに基づいて、推力発生装置１〜３の夫々のプロペラ推力及び舵２２，３２の夫々の舵角を演算する（ステップＳ４）。
【００６６】
そして、ＣＰＵ５１は、ステップＳ４にて演算した各プロペラ推力及び各舵角を、入出力インタフェース５４から制御装置９へ出力させ（ステップＳ５）、処理をステップＳ１へ戻す。
【００６７】
ステップＳ４の推力演算処理は、（２４）〜（２６）式のような推力のバランス関係、即ち船舶１０の動作モデルと、（３０）式のような評価関数とに基づいて定められた演算処理である。
【００６８】
【数４６】
2004042885.tif 000026
【数４７】
2004042885.tif 000027
【数４８】
2004042885.tif 000028
但し、ｕ１（ｔ）：推力発生装置１の正規化されたプロペラ推力
ｕ２（ｔ）：推力発生装置２の正規化されたプロペラ推力
ｕ３（ｔ）：推力発生装置３の正規化されたプロペラ推力
ｘ１（ｔ）：推力発生装置２の正規化された舵角
ｘ２（ｔ）：推力発生装置３の正規化された舵角
ｐ１１〜ｐ３５：所定の係数
η（ｘ）：（２７）式で与えられる関数
ν（ｘ）：（２８）式で与えられる関数
Ｘ：船舶１０の推力の前後方向成分
Ｙ：船舶１０の推力の左右方向成分
Ｎ：船舶１０の旋回方向のモーメント
【数４９】
2004042885.tif 000029
但し、ｋ１：定数
【数５０】
2004042885.tif 000030
但し、ｋ２：定数
【００６９】
また、推力発生装置２の舵２２の舵角速度ｕ４（ｔ）及び推力発生装置３の舵３２の舵角速度ｕ５（ｔ）は、（２９）式で与えられる。
【００７０】
【数５１】
2004042885.tif 000031
【数５２】
2004042885.tif 000032
但し、ｓｆ１，ｓｆ２：重み係数
ｒ１〜ｒ５：重み係数
ｑ１，ｑ２：重み係数
α：舵角の中立点を規定する定数
【００７１】
（２４）式は、次のようにして前述した（１）式から得られる。まず、推力発生装置２によって発生する推力のうち、船舶１０の前後方向の成分は、プロペラ２１によって発生する推力ｕ２（ｔ）から、舵２２によって生じる損失を差し引いたものとして概ね考えることができる。舵２２の舵角ｘ１（ｔ）が０の場合には、プロペラ２１によって後方へ送られた流体の流れに対し、舵２２は殆ど揚力を発生せず、プロペラ２１によって発生する推力ｕ２（ｔ）が、殆どそのまま推力発生装置２によって発生する推力の前後方向成分となる。一方、舵２２の舵角ｘ１（ｔ）が０以外の場合には、プロペラ２１によって後方へ送られた流体の流れに対し、舵２２に揚力が発生するとともに、その舵面に抗力が発生することとなる。従って、このとき、推力発生装置２の発生推力の前後方向成分からは、左右方向成分に対応する損失が生じる。
【００７２】
このような損失は、プロペラ２１が逆回転し、前方へ推力を発生した場合には生じない。このような推力発生装置２の非線形特性をモデル化するために、本実施の形態１では、前述した関数η（ｘ），ν（ｘ）を導入している。図４（ａ）は、関数η（ｘ）を示すグラフであり、図４（ｂ）は、関数ν（ｘ）を示すグラフである。関数η（ｘ）は、図４（ａ）に示すように、（３１）式で与えられる関数η１（ｘ）を、ｘ＝０の付近で滑らかになるように近似した関数である。
【００７３】
【数５３】
2004042885.tif 000033
【００７４】
また、関数ν（ｘ）は、図４（ｂ）に示すように、（３２）式で与えられる関数ν１（ｘ）を、ｘ＝０の付近で滑らかになるように近似した関数である。
【００７５】
【数５４】
2004042885.tif 000034
【００７６】
よって、（２４）式の左辺第１項は、ｕ２（ｔ）が負の場合には、ｐ１１ｕ２（ｔ）と略同値となり、ｕ２（ｔ）が正の場合には、ｐ１１［ｕ２（ｔ）−ｐｒ１ｕ２（ｔ）｜ｘ１（ｔ）｜］と略同値となる。ここで、ｘ１（ｔ）は、ｕ２（ｔ）ｘ１（ｔ）で、推力発生装置２による発生推力の左右方向成分（左右舷舵力）を表すような舵角とされているものとしている。また、定数ｐ１１は、プロペラ２１の形状等に基づいて定められる。
【００７７】
このように、（２４）式の左辺第１項は、推力発生装置２の発生推力の前後方向成分についての非線形特性をモデル化したものである。つまり、（２４）式は、（１）式においてｎを３とし、推力発生装置１の発生推力の前後方向成分を０としたものと考えることができ、船舶１０の推力の前後方向成分の釣り合い関係を、推力発生装置２，３の非線形特性に基づいてモデル化したものであるといえる。なお、推力発生装置３については、推力発生装置２と同様に考えることができるので、説明を省略している。
【００７８】
また、（２５）式は、（２）式のｎを３とし、ｆｙ１，ｆｙ２，ｆｙ３を（３３）〜（３５）式で与えたものである。
【００７９】
【数５５】
2004042885.tif 000035
【数５６】
2004042885.tif 000036
【数５７】
2004042885.tif 000037
【００８０】
前述したように、推力発生装置１は、プロペラが船舶１０の船体の横方向へ配置されており、舵等の揚力を発生する手段を具備していないため、その推力の殆どが船舶１０の左右方向に生じることとなる。従って、推力発生装置１のプロペラによって発生する推力ｕ１（ｔ）は、殆どそのまま推力発生装置１によって発生する推力の左右方向成分ｆｙ１となる。このため、推力発生装置１の発生推力の左右方向成分に関する特性は、（３３）式のように表すことができる。
【００８１】
一方、推力発生装置２には、舵２２が設けられているので、推力発生装置２の左右舷舵力ｕ２（ｔ）ｘ１（ｔ）に応じた推力が、推力の左右方向成分ｆｙ２となる。このため、推力発生装置２の発生推力の左右方向成分に関する特性は、（３４）式のように表すことができる。また、推力発生装置３も、推力発生装置２と略同様の構成であるので、その発生推力の左右方向成分に関する特性は、（３４）式と同様の（３５）式のように表すことができる。
【００８２】
また、（２６）式は、（３）式のｎを３とし、以下に説明するようなｋ２１，ｋ１２，ｋ２２，ｋ１３，ｋ２３が与えられたものである。まず、推力発生装置１に関して、（３）式及び（２６）式から、（３６）式が導かれる。
【００８３】
【数５８】
2004042885.tif 000038
【００８４】
そして、前述したような推力発生装置１の特性から、ｆｘ１＝０とすることができる。また、（３６）式と、（３３）式とによって、（３７）式のようなｋ２１が導かれる。
【００８５】
【数５９】
2004042885.tif 000039
【００８６】
次に、推力発生装置２に関して、（３）式及び（２６）式から、（３８）式が導かれる。
【００８７】
【数６０】
2004042885.tif 000040
【００８８】
また、（２４）式及び（３４）式より、次のような（３９）式が導かれる。
【００８９】
【数６１】
2004042885.tif 000041
【００９０】
従って、（３８）式及び（３９）式により、（４０）式及び（４１）式のようなｋ１２，ｋ２２が導かれる。
【００９１】
【数６２】
2004042885.tif 000042
【数６３】
2004042885.tif 000043
【００９２】
また、推力発生装置３の構成が、推力発生装置２の構成と略同様であることから、同様にして、（４２）式及び（４３）式のようなｋ１３，ｋ２３が導かれる。
【００９３】
【数６４】
2004042885.tif 000044
【数６５】
2004042885.tif 000045
【００９４】
なお、本実施の形態１においては、プロペラ２１及び舵２２によって推力発生装置２が構成され、プロペラ３１及び舵３２によって推力発生装置３が構成されたものとみなした場合についての（２４）〜（２６）式の求め方について説明したが、以下のように、プロペラ及び舵を夫々独立した推力発生装置とみなして、（２４）〜（２６）式を求めることもできる。
【００９５】
この場合では、プロペラ２１を推力発生装置２とし、プロペラ３１を推力発生装置３とし、舵２２を推力発生装置４とし、舵３２を推力発生装置５として考える。なお、舵はそれ自体で推力を発生することはないが、揚力等の舵が発生する作用力も一種の推力とみなすことができるため、ここでは舵を推力発生装置として扱うこととしている。
【００９６】
まず、各推力発生装置１〜５の発生推力の前後方向成分ｆｘ１〜ｆｘ５を、（４４）〜（４８）式で与える。
【００９７】
【数６６】
2004042885.tif 000046
【数６７】
2004042885.tif 000047
【数６８】
2004042885.tif 000048
【数６９】
2004042885.tif 000049
【数７０】
2004042885.tif 000050
【００９８】
ここで、推力発生装置４，５の発生推力の前後方向成分ｆｘ４，ｆｘ５は、前述した推力の損失として考えられるので、（４７）、（４８）式のように与えることができることがわかる。そして、（１）式と、以上の（４４）〜（４８）式とによって、（２４）式が求められる。
【００９９】
また、各推力発生装置１〜５の発生推力の左右方向成分ｆｙ１〜ｆｙ５を、（４９）〜（５３）式で与える。
【０１００】
【数７１】
2004042885.tif 000051
【数７２】
2004042885.tif 000052
【数７３】
2004042885.tif 000053
【数７４】
2004042885.tif 000054
【数７５】
2004042885.tif 000055
【０１０１】
ここで、推力発生装置４，５の発生推力の左右方向成分ｆｙ４，ｆｙ５は、前述した舵力として考えられるので、（５２）式及び（５３）式のように与えることができることがわかる。そして、（２）式と、以上の（４９）〜（５３）式とによって、（２５）式が求められる。
【０１０２】
図５は、推力発生装置１〜５の船舶１０上での位置関係を示す模式的平面図である。図５に示すように、船舶１０の船体中心を原点として、船首方向をζ方向、右舷方向をξ方向としたζ−ξ座標系を定義する。このとき、（ｋ２ｉ，ｋ１ｉ）（但し、ｉ＝１〜５）を、夫々推力発生装置ｉのζ−ξ座標系における座標値とする。そして、このようなｋ１１〜ｋ１５，ｋ２１〜ｋ２５，ｆｘ１〜ｆｘ５，ｆｙ１〜ｆｙ５と、（３）式とによって、（２６）式が求められる。
【０１０３】
このように、プロペラと舵とを別個の推力発生装置とみなすことによっても、（２４）〜（２６）式の推力のバランス関係を得ることができる。
ここで、（３０）式における舵角の中立点を規定する定数αについて説明する。図６は、舵２２，３２の状態を示す模式図である。船舶１０は、図６（ａ）に示すように、舵２２，３２を、略Ｖ字を形成するように夫々所定角度だけ内側へ傾斜させたり、若しくは図６（ｂ）に示すように、舵２２，３２を、略倒立Ｖ字を形成するように夫々所定角度だけ外側へ傾斜させることが可能である。以下、図６（ａ）の状態をＶ舵といい、図６（ｂ）の状態を逆ハ舵という。これらのように、舵２２，３２を夫々互いに反対方向へ同一角度だけ傾斜させた場合には、舵２２による左右舷舵力と、舵３２による左右舷舵力とが相殺され、船舶１０全体としては、推力発生装置２，３による左右方向への推力が発生されない。そこで、（３０）式では、Ｖ舵又は逆ハ舵の時の舵２２の舵角を正規化した値をαとすることとしている。これによって、Ｖ舵又は逆ハ舵の時の舵角の中立点を規定することができる。このようなＶ舵又は逆ハ舵を設定した場合には、プロペラ前進推力を発生させることにより、舵角速度が遅いときでも直ちに舵横力を発生させることができるという効果を奏する。
【０１０４】
更に、夫々のプロペラ推力ｕ１（ｔ），ｕ２（ｔ），ｕ３（ｔ）及び舵角速度ｕ４（ｔ），ｕ５（ｔ）には、次の（５４）式のような制限が設けられている。
【０１０５】
【数７６】
2004042885.tif 000056
但し、ｉ：１〜５の整数
ｕ１ｍ〜ｕ３ｍ：推力発生装置１〜３のプロペラ推力の上下限値
ｕ４ｍ，ｕ５ｍ：舵２２，３２の舵角速度の上下限値
【０１０６】
このように、推力発生装置１〜３のプロペラ推力、及び舵２２，３２の舵角速度には、所定の上下限値が存在するため、推力発生装置１〜３の動特性モデルは、（５４）式のような不等式の拘束条件を含むものとして表現される。このような不等式の拘束条件は、推力演算においてその取り扱いが困難であるため、（５５）式のような等式の拘束条件へ変換される。
【０１０７】
【数７７】
2004042885.tif 000057
但し、ｕｄ：ダミー変数
【０１０８】
また、（５５）式の等式拘束条件の代わりに、（５６）式のような等式拘束条件を採用することもできる。
【０１０９】
【数７８】
2004042885.tif 000058
【０１１０】
このような、ダミー変数を用いた等式拘束条件を採用する場合には、未知変数が増加し、この結果演算量の増大を招くこととなる。従って、単に演算によって得られたプロペラ推力又は舵角速度が上下限値を越えた場合には、強制的に上下限値をプロペラ推力又は舵角速度として設定するように構成することもできる。
【０１１１】
また、（３０）式の重み係数ｑ１，ｑ２の代わりに、（５７）式で表されるペナルティー法による重みｑｐ１，ｑｐ２を採用することもできる。
【０１１２】
【数７９】
2004042885.tif 000059
但し、ｉ：１又は２
ｑｉ：通常の重み係数
ｋｐ：係数
ｘｉｍ：舵２２，３２の舵角の上下限値
δ：正の定数
【０１１３】
なお、δは、δ＜ｘｉｍを満たすものとする。（５７）式の右辺第２項は、状態量ｘｉ（ｔ）を制限値ｘｉｍ−δ内に拘束するためのペナルティーであり、ｋｐはペナルティーの大きさを調節するためのパラメータである。これにより、ｑｐｉは、状態量ｘｉ（ｔ）が制限値ｘｉｍ−δ内にあるときには、通常の重み係数ｑｉと同値であり、状態量ｘｉ（ｔ）が制限値ｘｉｍ−δを越えたときにはｑｉに比べて非常に大きな値となる。また、ｍａｘ関数は、引数のうち大きい値を出力する関数であり、その二乗は微分可能な滑らかな関数となる。また、δは、重みｑｐｉを状態量ｘｉ（ｔ）がｘｉｍを越える前から大きくするように調整するためのパラメータである。
【０１１４】
また、ステップＳ２において、推力発生装置１〜３の故障の発生を検出した場合に、ステップＳ３でこの故障が発生した推力発生装置に係る重み係数の値を大幅に増加させることによって、この推力発生装置への推力発生の指令値が略０となり、他の正常な推力発生装置に対する推力配分を適切に行うことが可能となる。
【０１１５】
ステップＳ４の推力演算処理は、前述したような（２４）〜（２６）式の推力のバランス関係を満たしつつ、（３０）式の評価関数の値を最小化するプロペラ推力ｕ１（ｔ）〜ｕ３（ｔ），舵角速度ｕ４（ｔ），ｕ５（ｔ）及び舵角ｘ１（ｔ），ｘ２（ｔ）を求める処理である。
【０１１６】
ここで、ステップＳ４の推力演算処理の一例について説明する。なお、ここでは説明を簡単にするため、（２４）〜（３０）式を（１）〜（５）式及び（１１）式の形式で表現した場合について説明している。また、ｘ（ｔ），ｕ（ｔ）は、夫々（５８）式、（５９）式のようなベクトルである。また、（６０）式の関数Ｃを定義する。
【０１１７】
【数８０】
2004042885.tif 000060
【数８１】
2004042885.tif 000061
【数８２】
2004042885.tif 000062
【０１１８】
従って、（１）〜（３）式の等式拘束条件は、関数Ｃを用いて（６１）式のように表すことができる。
【０１１９】
【数８３】
2004042885.tif 000063
【０１２０】
まず、本例における推力演算処理では、評価区間ＴをＮステップに分割して離散近似された以下の問題を解くこととする。
【０１２１】
【数８４】
2004042885.tif 000064
【数８５】
2004042885.tif 000065
【数８６】
2004042885.tif 000066
【数８７】
2004042885.tif 000067
【０１２２】
ここで、ｘ＊ｊ（ｔ）は、ｘ（ｔ）を初期状態とする離散近似系のｊステップめの状態を表している。また、Δτ＝Ｔ／Ｎとしたときに、ｘ＊ｊ（ｔ）は、ｘ（ｔ＋ｊΔτ）である。
【０１２３】
このとき、（６６）〜（６８）式の条件を満たす必要がある。但し、下付き添字は偏微分を表すものとしている。
【０１２４】
【数８８】
2004042885.tif 000068
【数８９】
2004042885.tif 000069
【数９０】
2004042885.tif 000070
【０１２５】
但し、Ｈは（６９）式で定義されたハミルトニアンである。
【０１２６】
【数９１】
2004042885.tif 000071
【０１２７】
未知変数である制御入力の系列をまとめたベクトルＵを以下のように定義する。
【０１２８】
【数９２】
2004042885.tif 000072
【０１２９】
Ｕ（ｔ）とｘ（ｔ）とが与えられた場合、（６６）式は、（７１）式のようなｘ（ｔ）とＵ（ｔ）に関する方程式とみなせる。
【０１３０】
【数９３】
2004042885.tif 000073
【０１３１】
ここで、Ｇ（Ｕ，ｘ，ｔ）が０に収束するようにＵ（ｔ）を更新すべく、（７２）式を定義する。
【０１３２】
【数９４】
2004042885.tif 000074
但し、ζ：正の実数
【０１３３】
ヤコビアンＧＵが正則であれば、（７３）式の微分方程式を積分することにより、Ｕ（ｔ）を更新していくことができる。
【０１３４】
【数９５】
2004042885.tif 000075
【０１３５】
ここで、（７４）〜（７６）式で与えられるベクトルＷ，ｗ及び実数ωに対して、それらとヤコビアンとの積を（７７）式のように前進差分で近似する。
【０１３６】
【数９６】
2004042885.tif 000076
【数９７】
2004042885.tif 000077
【数９８】
2004042885.tif 000078
【数９９】
2004042885.tif 000079
但し、β：正の実数
【０１３７】
このとき、（７２）式は、次の（７８）式で近似され、これは（７９）式と等価である。
【０１３８】
【数１００】
2004042885.tif 000080
【数１０１】
2004042885.tif 000081
【０１３９】
但し、ｂは（８０）式で定義されたものとする。
【０１４０】
【数１０２】
2004042885.tif 000082
【０１４１】
従って、（７９）式に基づいてＵ（ｔ）を演算することが可能である。ここでは、その演算処理の一例について説明する。図７，８は、推力演算処理の処理手順の一例を示すフローチャートである。この演算処理について、図７，８を参照しながら説明する。まずＣＰＵ５１は、変数ｒ０に−ｂ（Ｕ，ｘ，ドットｘ，ｔ）を代入する（ステップＳ１０）。次に、ＣＰＵ５１は、変数ｖ１にｒ０／‖ｒ０‖を代入し（ステップＳ１１）、変数γにノルム‖ｒ０‖を代入する（ステップＳ１２）。また、ＣＰＵ５１は、変数εにγを代入し（ステップＳ１３）、カウンタｋを０にセットする（ステップＳ１４）。
【０１４２】
次にＣＰＵ５１は、ｋが、逆行列演算の繰り返し回数を表す所定の正の整数ｋｍａｘより小さいか否かを判別する（ステップＳ１５）。
ステップＳ１５において、ｋ＜ｋｍａｘを満たす場合、ＣＰＵ５１は、ｋを１だけインクリメントする（ステップＳ１６）。そして、ＣＰＵ５１は、ｖｋ＋１に、ＤβＧ（Ｕ，ｘ＋βドットｘ，ｔ＋β：ｖｋ，０，０）を代入する（ステップＳ１７）。
【０１４３】
次にＣＰＵ５１は、カウンタκを１にセットし（ステップＳ１８）、変数ｈκｋに、ｖｋ＋１Ｔｖκを代入する（ステップＳ１９）。また、ＣＰＵ５１は、ｖｋ＋１を、ｖｋ＋１−ｈκｋｖκに更新する（ステップＳ２０）。
【０１４４】
ＣＰＵ５１は、κ＜ｋを満たすか否かを判別し（ステップＳ２１）、κ＜ｋの場合には、κを１だけインクリメントして（ステップＳ２２）、ステップＳ１９へ処理を移す。
【０１４５】
また、ステップＳ２１において、κ≧ｋの場合には、ＣＰＵ５１は、変数ｈｋ＋１，ｋにノルム‖ｖｋ＋１‖を代入し（ステップＳ２３）、ｖｋ＋１を、ｖｋ＋１／‖ｖｋ＋１‖に更新する（ステップＳ２４）。
次に、ＣＰＵ５１は、（８１）式及び（８２）式で表されるｅ１及びＨｋに対して、‖εｅ１−Ｈｋｙｋ‖を最小にするｙｋを求める（ステップＳ２５）。但し、ｙｋは（８３）式で表されるベクトルである。
【０１４６】
【数１０３】
2004042885.tif 000083
【数１０４】
2004042885.tif 000084
但し、ｉ＞ｊ＋１のとき、ｈｉｊ＝０
【数１０５】
2004042885.tif 000085
【０１４７】
次に、ＣＰＵ５１は、γを‖εｅ１−Ｈｋｙｋ‖に更新し（ステップＳ２６）、ステップＳ１５へ処理を移す。ステップＳ１５において、ｋ≧ｋｍａｘの場合、ＣＰＵ５１は、ドットＵをＶｋｙｋとする（ステップＳ２７）。但し、Ｖｋは（８４）式で表されるベクトルである。
【０１４８】
【数１０６】
2004042885.tif 000086
【０１４９】
そして、ＣＰＵ５１は、ドットＵを実時間で積分して、Ｕを算出し（ステップＳ２８）、リターンする。
【０１５０】
本願発明者らは、以上で説明したような推力演算アルゴリズムを用いて推力配分のシミュレーションを行うコンピュータプログラムをコンピュータに実装し、推力配分のシミュレーションの実験を行った。図９（ａ）は、本実験で使用した各パラメータの設定値を示す図表であり、図９（ｂ）は、本実験で使用した各係数の設定値を示す図表である。図９（ａ）に示すように、本実験では、サンプリング周期を０．１秒とし、ζを１０．０とした。また、（２７）式におけるｋ１及び（２８）式におけるｋ２を夫々５０．０とした。また、評価区間Ｔを０．２秒とし、舵角の中立点を規定する定数αを０．５とした。更に逆行列演算の繰り返し回数ｋｍａｘを２０とし、評価区間Ｔの分割数Ｎを２とした。
【０１５１】
また、図９（ｂ）に示すように、本実験では、重み係数ｒ１〜ｒ５を、夫々０．１，１．０，１．０，０．１，０．１とし、重み係数ｓｆ１，ｓｆ１を夫々０．００１とした。また、重み係数ｑ１，ｑ２を夫々０．１とし、係数ｋｐを１０００、δを０．０２とした。
【０１５２】
また、本実験では、６０秒間だけ船舶１０の目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントを変化させた場合について、各推力発生装置１〜３のプロペラ推力ｕ１（ｔ）〜ｕ３（ｔ）並びに舵２２，３２の舵角ｘ１（ｔ），ｘ２（ｔ）及び舵角速度ｕ４（ｔ），ｕ５（ｔ）を演算した。また、演算した結果のｕ１（ｔ）〜ｕ５（ｔ）及びｘ１（ｔ），ｘ２（ｔ）から、船舶１０の全体の推力の前後方向成分及び左右方向成分並びに旋回方向のモーメントを演算し、夫々の目標値との比較を行った。更に、図７，８に示した推力演算処理で得られたｕ１（ｔ）〜ｕ５（ｔ）及びｘ１（ｔ），ｘ２（ｔ）が、オイラー・ラグランジュの条件式（１６）をどの程度満足しているかを示す指標として、条件式（１６）を満たしていない程、数値が大きくなる最適性の誤差を求めた。
【０１５３】
図１０，１１は、本実験の結果の一例を示すグラフである。図１０（ａ）は、推力発生装置１のプロペラ推力ｕ１（ｔ）の時間的変化を示すグラフであり、図１０（ｂ）は、推力発生装置２，３のプロペラ推力ｕ２（ｔ），ｕ３（ｔ）の時間的変化を示すグラフである。図１０（ｂ）において、実線はプロペラ推力ｕ２（ｔ）の変化を示しており、一点鎖線はプロペラ推力ｕ３（ｔ）の変化を示している。また、図１０（ｃ）は、舵２２，３２の舵角ｘ１（ｔ），ｘ２（ｔ）の時間的変化を示すグラフであり、図１０（ｄ）は、舵角速度ｕ４（ｔ），ｕ５（ｔ）の時間的変化を示すグラフである。図１０（ｃ）において、実線は舵角ｘ１（ｔ）の変化を、一点鎖線は舵角ｘ２（ｔ）の変化を示している。また、図１０（ｄ）において、実線は舵角速度ｕ４（ｔ）の変化を、一点鎖線は舵角速度ｕ５（ｔ）の変化を示している。
【０１５４】
図１１（ａ）は、船舶１０の全体の推力の前後方向成分の時間的変化を示すグラフである。また、図１１（ｂ）は、船舶１０の全体の推力の左右方向成分の時間的変化を示すグラフであり、図１１（ｃ）は、船舶１０の旋回方向のモーメントの時間的変化を示すグラフである。図１１（ｄ）は、推力演算処理での誤差の総和の時間的変化を示すグラフである。図１１（ａ）〜（ｃ）において、実線は、推力発生装置１〜３がｕ１（ｔ）〜ｕ５（ｔ）及びｘ１（ｔ），ｘ２（ｔ）で動作した場合の船舶１０の推力の前後方向成分及び左右方向成分並びに旋回方向のモーメントの時間的変化を夫々表しており、一点鎖線は、目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントの時間的変化を表している。
【０１５５】
図１１（ａ）〜（ｃ）に示すように、目標推力の前後方向成分及び左右方向成分並びに旋回方向の目標モーメントは、夫々異なる周期の正弦波状に変化するものとした。また、αを０．５としていることからわかるように、５０％のＶ舵とした場合について実験を行った。
【０１５６】
この結果、概ね目標値に一致する推力及びモーメントが得られた。また、推力演算処理での誤差の総和についても、十分に小さい誤差に抑えることができた。
【０１５７】
また、図１２，１３は、本実験の結果の他の例を示すグラフである。図１３（ａ）〜（ｃ）に示すように、本例においては、シミュレーション開始から１０秒後に、目標推力の前後方向成分及び左右方向成分を０．２として与え、３０秒後に推力発生装置１に故障が発生したものとしてシミュレーションを行った。また、故障発生を検出した後、推力発生装置１に係る重み係数ｒ１の値を、０．１から１０００へと変更した。
【０１５８】
この結果、推力発生装置１の故障発生時点たるシミュレーション開始３０秒後に、ｕ１（ｔ）が約０．２から０へと変化し、ｕ２（ｔ）が約０．２から約０．７へと変化し、ｕ３（ｔ）が約０．２から約−０．４へと変化した。また、舵角ｘ１（ｔ），ｘ２（ｔ）も、この時点から夫々若干変化した。
【０１５９】
そして、船舶１０の全体の推力の前後方向成分が、シミュレーション開始３０秒後から、目標値から若干離れているものの、全体としては概ね目標値に一致する推力及びモーメントが得られた。また、この場合の推力演算処理での誤差の総和は、シミュレーション開始３０秒後付近で一時的に大きくなったが、時間の経過とともに減少した。
【０１６０】
これにより、本発明に係る推力配分方法を用いることによって、推力発生装置に故障が発生した場合であっても、適切な推力配分を継続して行うことが可能であることがわかる。
【０１６１】
また、本実験では、標準的な構成のパーソナルコンピュータを用いたが、６０秒間に亘るｕ１（ｔ）〜ｕ５（ｔ）及びｘ１（ｔ），ｘ２（ｔ）を演算するのに、約４秒間しか要さなかった。このことから、本実施の形態１に係る推力演算方法を用いることにより、実時間での推力配分を十分に行うことが可能であることがわかる。
【０１６２】
（実施の形態２）
図１４は、本発明の実施の形態２に係る推力配分装置が搭載された船舶の構成の一例を示す模式図である。図１４に示すように、船舶１１０は、推力発生装置１０１〜１０６を備えている。かかる推力発生装置１０１〜１０６は、首振り式の推力発生装置、所謂旋回式スラスタであり、推力発生装置１０１〜１０６自体が旋回することによって、その推力発生方向を変化させることが可能である。
【０１６３】
なお、図１４においては、船舶１１０が６つの推力発生装置１０１〜１０６を船体の前部、中部、後部の左右両端に夫々有する構成としたが、推力発生装置の数及び配置位置は、これに限定されるものではない。また、以下においては、任意の数の旋回式の推力発生装置を船舶の任意の位置に配置した場合について説明する。
【０１６４】
図１５は、本発明の実施の形態２に係る推力配分装置が搭載された船舶１１０の構成を示すブロック図である。以下の説明においては、船舶１１０が、ｎ基の旋回式推力発生装置１０１〜１０ｎを有する場合について考える。図１５に示すように、船舶１１０は、ｎ基の推力発生装置１０１〜１０ｎを有しており、これらの他に、駆動装置１４１〜１４ｎ，推力配分装置１５０，検出装置１６０，設定装置１７０，目標推力演算装置１８０及び制御装置１９０を備えている。
【０１６５】
駆動装置１４１〜１４ｎは、夫々対応する推力発生装置１０１〜１０ｎに接続されている。これらの駆動装置１４１〜１４ｎは、例えばエンジン、変速機、及び油圧装置等から夫々構成されており、旋回式の推力発生装置１０１〜１０ｎのプロペラを、夫々所要の回転数によって回転させ、また推力発生装置１０１〜１０ｎ自体を、夫々所要の角度だけ旋回させるようになっている。
【０１６６】
検出装置１６０は、実施の形態１で説明した位置検出手段６１，速度検出手段６２及び方位角検出手段６３と同様の構成の位置検出手段１６１，速度検出手段１６２及び方位角検出手段１６３を備えている。
【０１６７】
また、設定装置１７０及び目標推力演算装置１８０の構成は、実施の形態１にて説明した設定装置７及び目標推力演算装置８の構成と同様であるので、その説明を省略する。
【０１６８】
推力配分装置１５０は、ＣＰＵ１５１，ＲＯＭ１５２，ＲＡＭ１５３，及び入出力インタフェース１５４から構成されており、入出力インタフェース１５４が目標推力演算装置１８０にデータを受信することが可能であるように接続されていて、目標推力演算装置１８０から入力した目標推力の前後方法成分及び左右方向成分並びに旋回方向の目標モーメントに基づいて、推力発生装置１０１〜１０ｎの夫々のプロペラ推力及び旋回角度を演算する。このような推力配分装置１５０は、入出力インタフェース１５４が制御装置１９０にデータを送信することが可能であるように接続されており、演算結果の各プロペラ推力及び各旋回角度を制御装置１９０に出力するようになっている。
【０１６９】
制御装置１９０は、駆動装置１４１〜１４ｎに接続されている。また、船舶１１０には、推力発生装置１０１〜１０ｎのプロペラの回転数を夫々検出する回転数センサと、旋回角度を夫々検出する旋回角度センサとが設けられている。これらの回転数センサ及び旋回角度センサは、制御装置１９０に接続されており、検出結果の回転数及び旋回角度の情報を、夫々制御装置１９０へ出力するようになっている。
【０１７０】
また、制御装置１９０は、ＣＰＵ、ＲＯＭ及びＲＡＭ等から構成されており、推力配分装置１５０から入力した各プロペラ推力及び旋回角度と、各種センサから入力した現在の各プロペラの回転数及び各旋回角度とに基づいて、旋回式推力発生装置１０１〜１０ｎの各プロペラ推力及び各旋回角度が、推力配分装置１５０から出力された各プロペラ推力及び各旋回角度に追従するように、例えばＰＩＤ制御等を用いて、推力発生装置１０１〜１０ｎの各プロペラ推力及び各旋回角度を制御する。この制御装置１９０から駆動装置１４１〜１４ｎに対して制御信号が出力され、推力発生装置１０１〜１０ｎが夫々駆動される。
【０１７１】
更に、制御装置１９０は、各種センサの測定値から得られた旋回式推力発生装置１０１〜１０ｎのプロペラ推力及び旋回角度の内の何れかが、推力配分装置１５０から入力したプロペラ推力及び旋回角度に追従しているか否かを判別し、追従していない場合には、その旋回式推力発生装置が故障したものとみなして、故障発生の情報を推力配分装置１５０へ出力するようにもなっている。
【０１７２】
次に、推力配分装置１５０が実行するプログラムの推力演算アルゴリズムの概要について説明する。図１６は、旋回式の推力発生装置の位置及び旋回角度について説明する模式図である。図１５に示すように、船体重心を原点とする船体固定水平座標系ｘｏ−ｙｏ（船首方向ｘｏ，右舷方向ｙｏ）に対して、ｉ番目の旋回式推力発生装置１０ｉの位置を（ｘｏｉ，ｙｏｉ）とし、プロペラ推力をＴｉ，ｘｏ軸からの旋回角度をθｉとすると、ｎ基の旋回式推力発生装置が発生する推力の船体前後方向成分Ｘ，左右方向成分Ｙ，旋回方向のモーメントＮは、（８）式、（９）式及び（１０）式で与えられる。
【０１７３】
【数１０７】
2004042885.tif 000087
【数１０８】
2004042885.tif 000088
【数１０９】
2004042885.tif 000089
【０１７４】
ここで、旋回式推力発生装置の状態変数ベクトルをｘ＝［Ｔ１，Ｔ２，…，Ｔｎ，θ１，θ２，…，θｎ］Ｔとし、制御装置１９０から各旋回式推力発生装置１０１〜１０ｎに与えられる操作量ｕをｕ＝ｄｘ／ｄｔとする。そして、前述した（１１）式にて表される評価関数Ｊを規定するスカラー値関数φ，Ｌを、夫々（８５）式及び（８６）式のように定義する。
【０１７５】
【数１１０】
2004042885.tif 000090
【数１１１】
2004042885.tif 000091
但し、ｘｒｊ：ｊ番目の状態変数の目標値
Ｘｒ：目標推力演算装置１８０から与えられた船舶１１０の目標推力の船体前後方向成分
Ｙｒ：目標推力演算装置１８０から与えられた船舶１１０の目標推力の船体左右方向成分
Ｎｒ：目標推力演算装置１８０から与えられた船舶１１０の旋回方向の目標モーメント
ｓｆｊ，ｓｗＸ，ｓｗＹ，ｓｗＮ，ｒｊ：重み係数
【０１７６】
また、操作量ｕ及び状態量ｘについて、（８７）式及び（８８）式のような拘束条件を定める。
【０１７７】
【数１１２】
2004042885.tif 000092
【数１１３】
2004042885.tif 000093
但し、ｕＬｊ：ｊ番目の操作量ｕｊの下限値
ｕＨｊ：ｊ番目の操作量ｕｊの上限値
ｘＬｊ：ｊ番目の状態量ｘｊの下限値
ｘＨｊ：ｊ番目の状態量ｘｊの上限値
【０１７８】
更に、旋回式の推力発生装置１０１〜１０ｎの何れかを故障等で停止させる場合には、対応する操作量の重み係数ｒｊを大きくするか、対応する操作量の下限値ｕＬｊ及び上限値ｕＨｊを０に近付けることで、当該推力発生装置の停止後の推力配分を正確に行うことが可能である。
【０１７９】
本願発明者らは、以上で説明したような推力演算アルゴリズムを用いて推力配分のシミュレーションを行うコンピュータプログラムをコンピュータに実装し、推力配分のシミュレーションの実験を行った。図１７は、実験１の概要を説明する模式図である。図１７に示す如く、本実験１においては、図１５で示した構成の船舶１１０に対して、右斜め前方から船体に向かう方向（船体に対して４５°の方向）の風速２５ｍ／ｓの風、及びこれと同一方向の潮力３ｋｔの潮流が外乱として与えられた場合に、船舶１１０に前方へ１０ｔの推力を発生させることを想定し、このときの１００秒間の各推力発生装置１０１〜１０６のプロペラ推力及び旋回角度の推移をシミュレートした。
【０１８０】
図１８は、本実験１の結果を示す模式的平面図であり、図１９〜図２５は、本実験１の結果を示すグラフである。また図１９は、推力発生装置１０１〜１０６のプロペラ推力の変化速度の時間的変化を示しており、図２０は、推力発生装置１０１〜１０６のプロペラ推力の時間的変化を示しており、図２１は、推力発生装置１０１〜１０６の旋回角速度の時間的変化を示しており、図２２は、推力発生装置１０１〜１０６の旋回角度の時間的変化を示しており、図２３は、船舶１１０全体の推力の前後方向成分の時間的変化を示しており、図２４は、船舶１１０全体の推力の左右方向成分の時間的変化を示しており、図２５は、船舶１１０全体の旋回推力（旋回方向のモーメント）の時間的変化を示している。
【０１８１】
図１８において、矢印は対応する箇所の推力発生装置１０１〜１０６の発生推力を示しており、矢印の向きでその発生推力の向き（即ち旋回角度）を示し、矢印の長さでその発生推力の大きさを示している。また、同じ箇所の矢印で上に重ねられて描かれているもの程、後で発生した推力を示すものとしている。即ち、図１８においては、時間の経過に従って、発生推力を示す矢印を上に重ねて描画することで、各推力発生装置１０１〜１０６の発生推力の向き及び大きさの時間的変化を示している。
【０１８２】
また、図１９〜２２において、推力発生装置１０１についての実験結果はピッチが長い破線で示しており、推力発生装置１０２についての実験結果は細い実線で示しており、推力発生装置１０３についての実験結果はピッチが短い破線で示しており、推力発生装置１０４についての実験結果は太い実線で示しており、推力発生装置１０５についての実験結果は一点鎖線で示しており、推力発生装置１０６についての実験結果は二点鎖線で示している。なお、本実験１では、各推力発生装置１０１〜１０６の最大推力を１７ｔとし、最小推力を最大推力の３０％としている。
【０１８３】
また、図２３〜２５においては、船舶１１０の目標推力の船体前後方向成分、船舶１１０の目標推力の船体左右方向成分、及び船舶１１０の旋回方向の目標モーメントを夫々一点鎖線にて示しており、実験結果の推力の船体前後方向成分、船体左右方向成分、及び旋回方向のモーメントを夫々実線にて示している。
【０１８４】
図１８〜２２に示すように、推力発生装置１０１〜１０６のプロペラ推力の変化速度、プロペラ推力、旋回角速度、及び旋回角度は、シミュレーション開始１０秒後から４０秒後までに亘って大きく変化し、その後は大きな変化がなく安定していることがわかる。
【０１８５】
また、図２３〜２５に示すように、船舶１１０全体の発生推力の前後方向成分及び左右方向成分並びに旋回方向のモーメントは、夫々シミュレーション開始１０秒後から４０秒後までに亘って漸増し、その後は略一定となっていることがわかる。更に、これは、夫々の目標値と略一致しており、船舶１１０全体の発生推力の前後方向成分及び左右方向成分並びに旋回方向のモーメントは、与えられた目標値にリアルタイムに追従していることがわかる。
【０１８６】
本願発明者らは、次のような実験２も行って、本実施の形態２で説明した推力演算アルゴリズムの有効性を確認した。図２６は、実験２の概要を説明する模式図である。図２６に示す如く、本実験２においては、図１５で示した構成の船舶１１０に対して、実験１と同様の条件で、船体左前側の推力発生装置１０１が故障により停止した場合を想定し、このときの１００秒間の各推力発生装置１０２〜１０６のプロペラ推力及び旋回角度の推移をシミュレートした。
【０１８７】
図２７は、本実験２の結果を示す模式的平面図であり、図２８〜図３４は、本実験２の結果を示すグラフである。また図２８は、推力発生装置１０２〜１０６のプロペラ推力の変化速度の時間的変化を示しており、図２９は、推力発生装置１０２〜１０６のプロペラ推力の時間的変化を示しており、図３０は、推力発生装置１０２〜１０６の旋回角速度の時間的変化を示しており、図３１は、推力発生装置１０２〜１０６の旋回角度の時間的変化を示しており、図３２は、船舶１１０全体の推力の前後方向成分の時間的変化を示しており、図３３は、船舶１１０全体の推力の左右方向成分の時間的変化を示しており、図３４は、船舶１１０全体の旋回推力（旋回方向のモーメント）の時間的変化を示している。
【０１８８】
図２８〜３１に示すように、推力発生装置１０２〜１０６のプロペラ推力の変化速度、プロペラ推力、旋回角速度、及び旋回角度は、シミュレーション開始１０秒後から４０秒後までの間は、実験１の結果よりも更に大きく変化し、その後は大きな変化がなく安定していることがわかる。また、シミュレーション開始４０秒後以降において、推力発生装置１０２〜１０６の発生推力が実験１の結果よりも概ね大きくなっており、推力発生装置１０１が停止している分が補われていることがわかる。
【０１８９】
また、図３２〜３４に示すように、船舶１１０全体の発生推力の前後方向成分及び左右方向成分並びに旋回方向のモーメントは、実験１の結果と略同様となっていることがわかる。このことから、実施の形態２に係る推力演算アルゴリズムを用いることによって、推力発生装置に故障が発生した場合であっても、適切な推力配分を継続して行うことが可能であることがわかる。
【０１９０】
【発明の効果】
本発明に係る推力配分方法及び推力配分装置によれば、操作量と発生推力との関係が線形性を示さない推力発生装置又は動特性を有する推力発生装置を航走体が備えている場合であっても、目的とする発生推力を得るために必要な操作量を正確に演算することができ、推力発生装置の発生推力を従来に比して高精度に制御することが可能となる。
【０１９１】
また、推力発生装置の故障等の不測の事態が発生した場合であっても、故障発生等に係る推力発生装置に関する重み係数をオンラインで変更することによって、目的とする航走体の推力を継続して得ることができる等、本発明は優れた効果を奏する。
【図面の簡単な説明】
【図１】本発明の実施の形態１に係る推力配分装置が搭載された船舶の構成を示す模式図である。
【図２】本発明の実施の形態１に係る推力配分装置が搭載された船舶の構成を示すブロック図である。
【図３】本発明の実施の形態１に係る推力配分装置の演算処理の流れを示すフローチャートである。
【図４】本発明の実施の形態１に係る推力発生装置の非線形特性のモデル化に用いられる関数を示すグラフであり、（ａ）は、関数η（ｘ）を示しており、（ｂ）は関数ν（ｘ）を示している。
【図５】プロペラ及び舵を別個の推力発生装置とみなした場合の各推力発生装置の船舶上での位置関係を示す模式的平面図である。
【図６】本発明の実施の形態１に係る船舶が備える舵の状態を示す模式図であり、（ａ）は、Ｖ舵を示しており、（ｂ）は、逆ハ舵を示している。
【図７】本発明の実施の形態１に係る推力配分装置の処理手順の一例を示すフローチャートである。
【図８】本発明の実施の形態１に係る推力配分装置の処理手順の一例を示すフローチャートである。
【図９】本発明の実施の形態１に係る推力配分装置について行った実験に用いたパラメータ及び係数の設定値を示す図表であり、（ａ）は、各パラメータの設定値を示しており、（ｂ）は、各係数の設定値を示している。
【図１０】本発明の実施の形態１に係る推力配分装置について行った実験の結果の一例を示すグラフであり、（ａ），（ｂ）は、各推力発生装置のプロペラ推力の時間的変化を示しており、（ｃ）は、各舵の舵角の時間的変化を示しており、（ｄ）は、各舵の舵角速度の時間的変化を示している。
【図１１】本発明の実施の形態１に係る推力配分装置について行った実験の結果の一例を示すグラフであり、（ａ）は、船舶全体の推力の前後方向成分の時間的変化を示しており、（ｂ）は、船舶全体の推力の左右方向成分の時間的変化を示しており、（ｃ）は、船舶全体の旋回方向のモーメントの時間的変化を示しており、（ｄ）は、本発明の実施の形態１に係る推力演算処理での誤差の総和の時間的変化を示している。
【図１２】本発明の実施の形態１に係る推力配分装置について行った実験の結果の他の例を示すグラフであり、（ａ），（ｂ）は、各推力発生装置のプロペラ推力の時間的変化を示しており、（ｃ）は、各舵の舵角の時間的変化を示しており、（ｄ）は、各舵の舵角速度の時間的変化を示している。
【図１３】本発明の実施の形態１に係る推力配分装置について行った実験の結果の他の例を示すグラフであり、（ａ）は、船舶全体の推力の前後方向成分の時間的変化を示しており、（ｂ）は、船舶全体の推力の左右方向成分の時間的変化を示しており、（ｃ）は、船舶全体の旋回方向のモーメントの時間的変化を示しており、（ｄ）は、本発明の実施の形態１に係る推力演算処理での誤差の総和の時間的変化を示している。
【図１４】本発明の実施の形態２に係る推力配分装置が搭載された船舶の構成の一例を示す模式的平面図である。
【図１５】本発明の実施の形態２に係る推力配分装置が搭載された船舶の構成を示すブロック図である。
【図１６】旋回式の推力発生装置の位置及び旋回角度について説明する模式図である。
【図１７】実験１の概要を説明する模式的平面図である。
【図１８】実験１の結果を示す模式的平面図である。
【図１９】実験１における推力発生装置のプロペラ推力の変化速度の時間的変化を示すグラフである。
【図２０】実験１における推力発生装置のプロペラ推力の時間的変化を示すグラフである。
【図２１】実験１における推力発生装置の旋回角速度の時間的変化を示すグラフである。
【図２２】実験１における推力発生装置の旋回角度の時間的変化を示すグラフである。
【図２３】実験１における船舶全体の推力の前後方向成分の時間的変化を示すグラフである。
【図２４】実験１における船舶全体の推力の左右方向成分の時間的変化を示すグラフである。
【図２５】実験１における船舶全体の旋回方向のモーメントの時間的変化を示すグラフである。
【図２６】実験２の概要を説明する模式的平面図である。
【図２７】実験２の結果を示す模式的平面図である。
【図２８】実験２における推力発生装置のプロペラ推力の変化速度の時間的変化を示すグラフである。
【図２９】実験２における推力発生装置のプロペラ推力の時間的変化を示すグラフである。
【図３０】実験２における推力発生装置の旋回角速度の時間的変化を示すグラフである。
【図３１】実験２における推力発生装置の旋回角度の時間的変化を示すグラフである。
【図３２】実験２における船舶全体の推力の前後方向成分の時間的変化を示すグラフである。
【図３３】実験２における船舶全体の推力の左右方向成分の時間的変化を示すグラフである。
【図３４】実験２における船舶全体の旋回方向のモーメントの時間的変化を示すグラフである。
【符号の説明】
１〜３　推力発生装置
２１，３１　プロペラ
２２，３２　舵
４１，４２　駆動装置
５　推力配分装置
５１　ＣＰＵ
５２　ＲＯＭ
５３　ＲＡＭ
５４　入出力インタフェース
６　検出装置
６１　位置検出手段
６２　速度検出手段
６３　方位角検出手段
７　設定装置
８　目標推力演算装置
９　制御装置
１０　船舶
１０１〜１０ｎ　推力発生装置
１４１〜１４ｎ　駆動装置
１５０　推力配分装置
１６０　検出装置
１６１　位置検出手段
１６２　速度検出手段
１６３　方位角検出手段
１７０　設定装置
１８０　目標推力演算装置
１９０　制御装置
１１０　船舶
【図１】
2004042885.tif 000116
【図２】
2004042885.tif 000117
【図３】
2004042885.tif 000118
【図４】
2004042885.tif 000119
【図５】
2004042885.tif 000120
【図６】
2004042885.tif 000121
【図７】
2004042885.tif 000122
【図８】
2004042885.tif 000123
【図９】
2004042885.tif 000124
【図１０】
2004042885.tif 000125
【図１１】
2004042885.tif 000126
【図１２】
2004042885.tif 000127
【図１３】
2004042885.tif 000128
【図１４】
2004042885.tif 000129
【図１５】
2004042885.tif 000130
【図１６】
2004042885.tif 000131
【図１７】
2004042885.tif 000132
【図１８】
2004042885.tif 000133
【図１９】
2004042885.tif 000134
【図２０】
2004042885.tif 000135
【図２１】
2004042885.tif 000136
【図２２】
2004042885.tif 000137
【図２３】
2004042885.tif 000138
【図２４】
2004042885.tif 000139
【図２５】
2004042885.tif 000140
【図２６】
2004042885.tif 000141
【図２７】
2004042885.tif 000142
【図２８】
2004042885.tif 000143
【図２９】
2004042885.tif 000144
【図３０】
2004042885.tif 000145
【図３１】
2004042885.tif 000146
【図３２】
2004042885.tif 000147
【図３３】
2004042885.tif 000148
【図３４】
2004042885.tif 000149