計測自動制御学会東北支部 第 332 回研究集会 (2021.5.26)
資料番号 332-3多項式モータトルクモデルにもとづいたオブザーバによるプロペラトルク推定
Observer-based propeller torque estimation based on polynomial motor
torque model○阿部直人*，小野浩平**，佐藤淳***○ Naoto Abe，Kohei Ono，Atsushi Satoh*岩手大学大学院，**（株）グレープシステム, ***岩手大学*Iwate University Graduate School，**GRAPE SYSTEMS INC. ***Iwate University
キーワード : 双線型オブザーバ，プロペラトルク，Carleman linearization
連絡先 : 〒 020-8551 岩手県盛岡市上田 4-3-5
岩手大学大学院 総合科学研究科 理工学専攻 機械・航空宇宙コース
佐藤1.

淳，Tel: 019-621-6404，E-mail: satsushi@iwate-u.ac.jp

はじめにが保障されない．線型システムを拡張した概念

近年，無人航空機 (UAV:Unmanned Air Vehi-cle) の利用が活発化している．機体としては，複数の独立したロータで発生させた推力やプロペラトルクの差を操舵力に利用するマルチコプタが主流である．
マルチコプタの飛行制御には，機体のダイナミクスと各プロペラ駆動系で発生する推力やプロペラトルクを知ることが重要である．プロペラのトルク特性はプロペラの回転速度 ω の２乗に比例する非線型な関係である 7) ．また，軸受けの摩擦特性もプロペラの回転速度に対して非線型であるため，プロペラの回転運動は非線型なダイナミクスをもつ．
非線型システムの解析や設計の難しさを避けるために，Taylor 展開に基づく線型化がしばしばおこなわれる．しかし，線型近似システムは，に双線型システムがある．双線型システムとは，状態空間表現が入力と状態の双線型項を含むシステムである．双線型制御理論は線型制御理論との共通部分が多い 2) ため，双線型システムは広範な応用が見込まれている．また，非線型システムの一種である多項式システムは，Carlemanlinearization により有限次元で双線型システムへ近似することができる．
本研究では，Carleman linearization を利用したオブザーバに基づき，プロペラトルクの推定をおこなう．プロペラトルクはプロペラ回転速度の多項式で表されるため，Carleman linearizationを利用した線型化が応用しやすい．オブザーバ
構成のために，プロペラ回転運動のモデリングと，運動モデルの Carleman linearization に基づいて有限次元近似をおこなう．

線型化する際に選んだ動作点近傍でしか有効性

–1–

2.

準備

2.1 問題設定
今回，駆動系は図 1 に示される岩手大学 佐藤研究室で開発中のクアッドコプタ MACTech に搭載されているプロペラ駆動系を想定する．よって MACTech に使用されているプロペラ，ブラシレス DC モータ，ESC を想定する．

Fig. 2

Fig. 1マルチコプタ推進系

クアッドコプタ MACTech

本研究で想定するマルチコプタの推進系を図2 に示す．推進系の運転状態を無風時のホバリング状態とし，プロペラの回転速度の範囲を 120 <=Fig. 3 ブラシレス DC モータのトルク‐スロットル‐回転数特性

ω<= 590[rad/s] と仮定する．
図 2 に示すように，モータのシャフトにはモー2.2タの出力である軸トルク Qs ，プロペラによる反トルクであるプロペラトルク Qp が作用する．軸
トルク Qs はモータのスロットルに応じてモータ内部で生ずるトルク成分値 Qm と，モータの軸受け摩擦による損失トルクである摩擦トルク

行列 A を m × n，行列 B を p × q の定数行列とする．A と B の Kronecker 積を以下の様に定義する 3) ．

[ 定義 1 ]成分 Qf からなるとする．すなわちQs (tn , ω) = Qm (tn ) + Qf (ω)

Kroneker 積(1)である．モータトルク Qm は tn スロットルの値 tn ，摩擦トルク Qf はプロペラの回転速度 ωに非線型に依存する．
そのため，Qs を tn と ω に関する多項式としてモデリングする．実験によりモデリングに利用するモータトルク特性は計測されて，そのデータを用いて MATLAB Curve ﬁtting toolbox で係数のフィッティングを行う．

a1n B
a2n B .. . 
am1 B am2 B . . . amn B

a11 B
 a21 BA⊗B =
 ..
 .

a12 B
a22 B..
.

...
...
..
.

(2)

このとき，A ⊗ B は mp × nq 行列となる．
また，x を n 次元状態ベクトルとすると，x[i]は次に定義するような x の i 次形式ベクトルである．

[ 定義 2 ]x = [x1 , x2 , . . . , xn ]T を n 次元ベクトル，p を–2–自然数とする．このとき，状態ベクトル同士のの時間微分が (6) 式で求められるとき，システムKronecker 積を以下のように定義する 3) ．

(5) は無限次元線型方程式に変換される．これをもちいた (5) 式の線型化を Carleman linearizationx ⊗ x = [x1 x1 , x1 x2 , . . . , x1 xn ,と呼ぶ．

. . . x n x1 , . . . x n xn ]T(3)ż = M z(8)

また，状態ベクトル同士の Kroneker 積により
ここで，定数行列 M の構造は (6) 式で計算さ生成されたベクトルの要素が，x の要素の p 次れた Aij をもちいて単項式:

xp11 xp22 . . . xpnn ,p=n∑pii=1を重複することなく全て辞書式順序で並べた

Npn 次元ベクトルを x[p] と書き，次式で定義したものを x の p 次形式と呼ぶ 4) ．ただし，Npn =n+p−1 Cp であり，p1 , p2 , . . . , pn は自然数または
A10 A11 . . . A1m
00 ...

0 . . . 0 A20 . . . A2m−1 A2m

M =
0 . . . A3m−2 A3m−1 A3m . . .
0

..
..
..
.
..
..
..
.
.
.
.
.
.
.
.
で表される．また，システム (5) に制御入力

u = [u1 , u2 , . . . ur ]T が加わる時，システム (5) は0 であり，x[0] := 1 と定義する．
x[p] := x[p−1] ⊗ xẋ = A0 + A1 x + A2 x[2] + . . . + An x[m] + Bu
(9)(4)となる．システム (9) の Carleman linearization

2.3はCarleman linearization

次式の状態方程式で表される n 次元 m 次多項ż = M0 z +式システムを考える．

1, …, n とする．いま，x の i 次形式ベクトルの時間微分がAij x[j+i−1](6)

j=0で表されるとき，(5) 式の Carleman linearizationを次のように定義する 3) ．ここで，Aij は In をn 次元の単位行列とするとき
Aij = A1j ⊗ In[i−1] + In ⊗ Ai−1j(7)で求められる定数行列である．

[ 定義 3 ][(10)のように双線型システムとなる 5) ．ここで各(5)ここで x ∈ Rn は状態ベクトル，Ai ∈ Rn×n , i =m∑

Mi zui + N0 ui=1ẋ = A0 + A1 x + A2 x[2] + . . . + An x[m]ẋ[i] =r∑

係数行列は (7) 式をもちいて以下のように計算される 6) ．

A10 A11 . . . A1m
00 ...

0 . . . 0 A20 . . . A2m−1 A2m

M0 = 
333

00.
.
.
A
A
A.
.
.
m
m−2m−1


..
..
..
.
.
.
..
..
..
..
.
.
.
.
 i,1
A1 Ai,1. . . Ai,100 ...
m2
 i,2i,2i,2i,20 . . .
A0 A1 . . . Am−1 AmMi = i,3i,3 0 Ai,3. . . Ai,3. . .0m−2 Am−1 Am

..
..
..
..
..
..
..
.
.
.
.
.
.
.
 B
 0
N =0 ..
.

n 次元状態ベクトルから生成した無限列 z =]Tx, x[2] , . . . に対して，システム (5) の i 形式で–3–

3.

プロペラトルク推定構造

3.3 システムの Carleman linearization

3.1 プロペラ回転運動のモデリング

変数を x1 := ω ，x2 := Qp とすると，(13) 式プロペラの回転運動のダイナミクスは以下のと (14) より，ようにあたえられる．

Jr ω̇ = Qs (tn , ω) − Qp

Jr ẋ1 = Qm0 + (k2 + k4 tn + k7 t2n )x1+ (k5 + k8 tn )x21 + k9 x31 − x2(11)1ẋ2 = − x2 + uvτ

ここで Jr は回転系の慣性モーメントである．
軸トルク Qs は，(15)となり，システム (15) は状態ベクトルが 2 次元

Qs (tn , ω) = k1 tn + k2 ω + k3 t2n + k4 tn ω + k5 ω 2+ k6 t3n + k7 t2n ω + k8 tn ω 2 + k9 ω 3 (12)と置く．(12) 式を (11) 式に代入して，起動トルでシステムの最大次数が 3 の多項式非線形システムである．

(15) 式を整理して，昇べきの順に並べなおすと，

ク Qm0 = k1 tn + k3 t2n + k6 t3n として整理すると

ẋ1 =Jr ω̇ = Qm0 + (k2 + k4 tn + k7 t2n )ω

k2 + k4 tn + k7 t2n
1x1 − x2JrJr
k5 + k8 tn 2 k9 3 Qm0x1 + x1 ++JrJr
Jr+ (k5 + k8 tn )ω 2 + k9 ω 3 − Qp (13)1ẋ2 = − x2 + uv
τとなり，これがプロペラの回転運動のダイナミクスを多項式システムでモデリングしたもので(16)となる．ここで，各係数をある．

k2 + k4 tn + k7 t2n1
, η2 =,
Jr2Jr
k 5 + k 8 tnk9
η3 =, η4 = ,JrJr
11
η5 = , η 6 = −Jrτη1 =3.2 プロペラトルのダイナミクスオブザーバをもちいてプロペラトルク Qp を推定するために，仮想的なダイナミクスQp+ uv
Q̇p = −τとし，入力の項を u1 = Qm0 ，u2 = uv とする．

(14)

また，測定可能な状態量は x1 = ω のみである
として，出力方程式と組み合わせると (16) 式はを導入する．

ẋ1 = η1 x1 + η2 x2 + η3 x21 + η4 x31 + η5 u1

本来プロペラトルク Qp の時間変化は様々な運転状態や環境要因等により非線形性を持つと

ẋ2 = η6 x2 + u2考えられるが，本研究では (14) 式のように Qpy = x1を線型信号モデル 7) であるとしてモデリングをおこなった．ここで，τ は Qp が取りうる値の帯域を制限するための正の時定数，uv は実際のプロペラトルクとのモデル化誤差をあたえる仮想(17)となる．

(17) 式の Carleman linearization をおこなうために，新たな状態変数としてz = [z1 , z2 , z3 , z4 , z5 , z6 , z8 , z9 ]T[]T= x1 , x2 , x21 , x1 x2 , x22 , x31 , x21 x2 , x1 x22 , x32的な入力である．

–4–

を定義する．z の要素の時間微分を計算し，推変数ベクトルは u = [u1 , u2 ]T である．



定対象の最高次数に合わせて 4 次以上の項を打ち切るとż1 = ẋ1 = η1 x1 + η2 x2 + η3 x21 + η4 x31 + η5 u1= η 1 z1 + η 2 z2 + η 3 z3 + η 4 z6 + η 5 u 1ż2 = ẋ2 = η6 x2 + u2 = η6 z2 + u2ż3 = 2x1 ẋ2≃ 2η1 x21 + 2η2 x1 x2 + 2η3 x21 + 2η5 x1 u1= 2η1 z3 + 2η2 z4 + 2η3 z6 + 2η5 z1 u1
ż4 = x1 ẋ2 + ẋ1 x2≃ (η1 + η6 ) x1 x2 + η2 x22 + η3 x21 x2 + η5 x2 u1 + x1 u2= (η1 + η6 ) z4 + η2 z5 + η3 z7 + η5 z2 u1 + z1 u2(18)となり，以降同様に z の時間微分を計算することで，推定対象 (17) 式は次の状態空間表現を得る．

ż = A0 z +2∑(20)
Ai zui + Bu,u = [u1 , u2 ]T

i=1

y = Czη1 η2 η300
η4
000
 0 η6 0
 0 0 2η12η20 2η3
0 00 η1 + η6 η20
A0 =  0 0002η6 0
0 0
000 3η1

0 0
0000


0000
0 0
0 0
0000

000

000 

000 

η300 
000 
3η200 2η1 + η62η20 η1 + 2η6 η2 0003η6(19)

ただし，各係数行列は以下のようになり，入力

000000000
 0
2η5 0000 0 η0005
A1 =  00000
 00 3η5 00
 000 2η5 0
000 η5
 000000000000000000000000000000000

0
0
0
0

0
0
0

0
0
(21)

–5–

また，(36) 式で求められる変換行列 T を用い

0001
A2 = 00000[B=000020000000000100000000020000000003000000000000000000000000000

η5 0 0 0 0 0 0 0 00 1 0 0 0 0 0 0 0

て新たな状態ベクトル ζ = T z を定義する．こ00
0
000000の ζ を用いて (24) 式を

x = Gζ + Hym(25)のように表現する．ここで，G，H は適当な定数行列である．これらの関係より次の定理が導ける．

[ 定理 1 ]システム (22) に対して，ym と u を入力とする]T

新たなシステム∑˙ζ̂ = F ζ̂ +T Ai V1 ym ui + T Bu + Lymri=1



1 0 0 0 0 0 0 0 0
C = 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0x̂ = Gζ̂ + Hym
ym = [x1 , x21 , x31 ]T(26)を考える．各係数行列 F ，G，H ，T ，L が3.4 双線型オブザーバの構成
Joachim の手法 8) をもちいてシステムの双線型オブザーバの設計をおこなう．状態変数ベクż = A0 z +(27)

T A − F T = LC(28)

T Ai V2 = 0, i = 1, 2

トル z ∈ R9 で近似された双線型システムは2∑Re (λ (F )) < 0GT + HC = E
Ai zui + Bu(30)を満足するとき，システム (25) は

i=1y = Cz(29)(22)

lim |x − x̂| = 0t→∞である．
測定可能な状態を ym ，測定不可能な状態をyu とすると，z は行列 V1 ，V2 をもちいて[ ][] ym
= V 1 y m + V2 yu(23)z = V 1 V2yu(31)となるような推定誤差 e = T z − ζ が漸近安定となるような，線型な誤差ダイナミクスė = F e(32)をもつ．ここで，λ (F ) は行列 F の固有値であと表せる．また，z = [x, x[2] , x[3] ]T , x = [x1 , x2 ]Tる．この係数行列 F ，G，H ，T ，L は次の定理であるため，変換行列 E をもちいてにより求められる．

[ 定理 2 ]x = Ezn × n 行列 F が行列 A0 とは異なる個別の固有値λ をもち，F の左固有ベクトル W をもちいて[E = I2]2×90(24)–6–( )F = W −1 diag λ W(33)で表されると仮定する．さらに対 (C, A) が可観である．また，変換行列 T も (36) 式よりパラ測であるとする．このとき，メータベクトルをもちいて計算される．これに T
p1−1  .. 
L=W  . より，

T
2.565 × 10−2
 6.028 × 10−1 −5.463 × 10−9 
−2.639 × 10−7 
T = −6.544 × 10−6 −2.235 × 10−7 
 2.658 × 10−14  1.357 × 10−12 
3.564 × 10−11(34)pTnに対する Sylvester 方程式

T A0 − F T = LC(35)の解 T は)−1 (
pT1 C A0 − λ1 In..

T = W −1 .


()
−1
T
p n C A0 − λ n I n(36)

pTµ ，µ = 1, . . . n は自由に割り当てられる h 次元パラメータベクトルである．このパラメータベクトルは次の行列方程式を解くことにより得られる．

()
pTµ C A − λµ Ai V2 = 0T ,i = 1, . . . , rとなる．最後に (41) 式により求めた T を (38) 式に代入し，一般化逆行列をもちいて G，H を求であたえられる．ここで，rank (T ) = n であり，µ = 1, . . . , n(41)(37)

また，G，H は T をもちいた線型方程式[ ][] T=E(38)
G HCめ．ここで，本研究における状態ベクトルの次元は 2 次元で，Carleman linearization 後の次元は 9 であるので，E は][1 0 0 0 0 0 0 0 0
E=0 1 0 0 0 0 0 0 0(42)となる．したがって，G，H は[]5.111 × 10−21
G=1.659
T1−4.256 × 10−2H = −5.111 × 10−22 9.063 × 10−9 −5.277 × 10−21 3.707 × 10−7(43)となる．

により求められる．

したがって，システム (17) のオブザーバ (26)

オブザーバ固有値を λ = −40 とすると，(33)式より F = −40，W = 1 となる．また，パラメータベクトルは (37) 式により求められる．本研究では[]
p = 1 −2.074 × 10−7 −8.169 × 10−6(39)を採用した．オブザーバゲイン L は (34) より．

W = 1 では L = pT である．したがって，[]
L = 1 −2.074 × 10−7 −8.169 × 10−6(40)はT
5.575 × 10−5˙ζ̂ = −40ζ̂ + −3.421 × 10−3  ym u10T[]T2.638 × 10−72
1.309×10+ 2.658 × 10−14  ym u2 +u6.018 × 10−10[]
5.111 × 10−17x̂ =ζ̂1.659
T
1−4.256 × 10−2+ −5.111 × 10−22 9.063 × 10−9  ym−5.277 × 10−21 3.707 × 10−7–7–



となる．

and Experiments”, IEEE Transactions on Control Systems Technology, Volume: 17, Issue:

4.

1, Jan. 2009

おわりに本研究では，マルチコプタの推進系を想定し，モータトルク特性の多項式モデルを用いて，プロペラトルクを推定するための双線型オブザー

参考文献
1) L.Pivano, et. al.，”A Four-Quadrant Thrust
Estimation Scheme for Marine Propellers:
Theory and Experiments”, IEEE Trans. on
Cont. Sys. Tech., Vol. 17, Issue: 1, Jan.
(2009)，1-12
2) Ronald R. Mohler，“Bilinear Control Processes”，ACADEMIC PRESS，(1973)，1-11
Kowalski，Willi-Hans

Steeb，“Nonlinear Dynamical Systems and Carleman
Linearization”，World

input-output linearization using Carleman linearization”，IEEE，2003 European ControlConference，(2003)

バを Carleman linearization により構成した．

3) Krzysztof

8) Joachim Deutscher，“Asymptotically exact

Scientiﬁc，(1991)，

7-8，11-12，83-93
4) 加藤誠，嘉納秀明，増淵正美，「多項式非線形系に対する状態観測器」，計測自動制御学会論文集，Vol.17，No.1，(1981)，8-15

5) J.Krener，“Linearization and bilinearizationof control systems”，Proceedings of the Allerton conference on circuit and system theory，(1974)
6) Negar Hashemian，Antonios Armaou，“FastMoving Horizon Estimation of Nonlinear Processes via Carleman Linearization”，American Control Conference，(2015)，3379-3385
7) Luca Pivano, Tor Arne Johansen and Øyvind
N. Smogeli ”A Four-Quadrant Thrust Estimation Scheme for Marine Propellers: Theory–8–