2 乗和多項式に基づくセミアクティブサスペンションの非線形制御
M2015SC010 松岡浩司
指導教員：大石泰章

1

はじめに自動車のセミアクティブサスペンションは，バネ上 (車体) を支えて衝撃を吸収するバネとバネの振動を減衰する減衰力可変式ショックアブソーバから構成され，減衰力可変式ショックアブソーバの減衰係数を調整することによりバネ上 (車体) の振動を抑制する．しかし，セミアクティブサスペンションは平衡点で入力が働かない双線形システムとなるため，制御することは難しい．双線形性を考慮しない制御系の設計法の 1 つとして，プラントを
図 1 セミアクティブサスペンション
アクティブサスペンションと仮定して線形制御理論を適用する方法があるが，乗員に制御の違和感を与えてしま
表 1 パラメータう．一方，双線形性を考慮する制御系の設計法の 1 つと
記号
詳細値して，パラメータ依存線形行列不等式 (PDLMI) を制約
Mb
バネ上 (車体) 質量375[kg]条件にもつロバスト半正定値計画問題 (SDP) を解く方法
Mw
バネ下 (車輪) 質量38[kg]がある．ロバスト SDP の解法の 1 つとして，PDLMI をKs
サスペンションのバネ定数
21560[N/m]
有限個の線形行列不等式 (LMI) へ帰着する方法 [1] が広Kt
タイヤのバネ定数
240000[N/m]
く知られているが，この場合の PDLMI はパラメータにCs
線形減衰係数
2000[Ns/m]
関して 2 次以上であるため，保守性が高くなってしまう．
Cv
可変減衰係数 −500〜500[Ns/m]
近年，ロバスト SDP の保守性の低い解法の 1 つとして，xb
バネ上 (車体) 変位[m]
2 乗和多項式 [2][3] に基づく手法が研究されている．
xw
バネ下 (車輪) 変位[m]本論文の構成は以下のとおりである．2 章では，セミアxr
路面変位[m]クティブサスペンションの説明をする．3 章では，擬似線形表現 [4] を導出する．また，コスト保証制御問題 [5] 3 制御系の設計を PDLMI を制約条件にもつロバスト SDP へ帰着させ， 3.1 擬似線形表現の導出
2 乗和多項式に基づきこのロバスト SDP の近似解を求め運動方程式 (1) 式と (2) 式に対して，状態変数を x(t) =る．4 章では，制御系の設計結果を示す．5 章では，数値 [x − x x − x ẋ ẋ ]T ，外乱入力を w(t) = ẋ ，制御rw wb w br
シミュレーションにより適用した方法の有用性を示す．6 入力を u(t) = C ，ストローク速度 ẋ − ẋ を θ と書く
vwb
章では，本論文のまとめと今後の課題を述べる．と，このシステムの擬似線形表現は以下のように求まる：

2{ẋ(t) = Ax(t) + B1 w(t) + B2 (θ)u(t)，

制御対象

本論文で扱うセミアクティブサスペンションを図 1 に示す．また，各パラメータは文献 [6] に従うものとする(表 1)．図 1 より，車両のバネ上 (車体)Mb とバネ下 (車輪)Mw の間にバネ Ks と減衰力可変式ショックアブソー
バ Cs + Cv からなるサスペンションが接続されている．また，路面変位 xr はタイヤ Kt を通してバネ下 (車輪)Mwに作用している．このモデルの運動方程式は，

Mb ẍb =Cs (ẋw − ẋb ) + Cv (ẋw − ẋb )+Ks (xw − xb )，(1)

Mw ẍw =−Cs (ẋw − ẋb ) − Cv (ẋw − ẋb )−Ks (xw − xb ) + Kt (xr − xw )(2)である．ここで，制御入力は可変減衰係数 Cv であり，状態変数の 1 つであるストローク速度 ẋw − ẋb と掛け算されている．そのため，ẋw − ẋb = 0 のときに入力が働かない双線形システムとなることが分かる．

y(t) = Cx(t) + D(θ)u(t)．
ここで，[]A B1 B2 (θ)
00−1 001=  KtKs
Cs
 Mw − Mw − Mw
Ks
Cs0Mb
Mb
[] [KsC D(θ) = 0 M
b0−1Cs
Mw
Cs−Mb
Cs
Mb

(3)

1000 
，0 − Mθw 
θ
0Mb]

Cs−Mb

θ
Mb

であり，B2 (θ) と D(θ) は θ に関して１次式である．また，
θ は −0.5 ≤ θ ≤ 0.5，θ の速度 θ̇ は −70 ≤ θ̇ ≤ 70 と仮定し，有界閉集合

Θ = {(θ, θ̇) ∈ R2 | − 0.5 ≤ θ ≤ 0.5, −70 ≤ θ̇ ≤ 70} (4)を定める．

3.2

ロバスト半正定値計画問題

3.4

擬似線形表現 (3) 式に対して，コスト保証制御問題の有解条件∫ ∞J=x(t)T Qx(t) + u(t)T Ru(t)dt ≤ γxT
0 x0 (5)0を達成する，状態フィードバック形式のコントローラ

u(t) = K(θ)x(t)(6)を設計する問題を考える．ここで，Q = QT = QT
h Qh ≻
O，R = RT ≻ O，γ > 0 であり，≻ (⪰) は対称行列が正定値 (半正定値) であることを表す．この問題の可解条件
は，PDLMI を制約条件にもつ以下のロバスト SDP を満足する決定変数 X(θ) = X(θ)T ⪰ O，F (θ)，γ が存在することである：

2 乗和多項式に基づく近似解法

領域 Θ が (θ, θ̇) の多項式不等式で規定される有界閉集合である場合，2 乗和多項式に基づきロバスト SDP(7)式を満たす近似解を求めることができる．準備として，有界閉集合 (4) 式を多項式 r1 (θ) = (θ + 0.5)(0.5 − θ)，
r2 (θ̇) = (θ̇ + 70)(70 − θ̇) により，

Θ = {(θ, θ̇) ∈ R2 |r1 (θ) ≥ 0, r2 (θ̇) ≥ 0}と書く．このとき，ロバスト SDP(7) 式の PDLMI 制約条件の十分条件として，以下を満たす多項式行列 S1 (θ)，S2 (θ)，S3 (θ̇)，S4 (θ) が存在する：

M1 (ξ, θ) − r1 (θ)S1 (θ) は 2 乗和多項式，
M2 (ξ, θ, θ̇) − r1 (θ)S2 (θ) − r2 (θ̇)S3 (θ̇) は 2 乗和多項式，最小化 γ

M3 (ξ, θ) − r1 (θ)S4 (θ) は 2 乗和多項式，条件

S1 (θ) は 2 乗和多項式，

M1 (ξ, θ) := X(θ) ⪰ O，−(Φ(θ) + Φ(θ)T − Ẋ(θ)) ∗
M2 (ξ, θ, θ̇) := Qh X(θ)
IRF (θ)0[]γI∗
M3 (ξ, θ) :=⪰O
I X(θ)(10)

S2 (θ) は 2 乗和多項式，∗∗  ⪰ O，R

S3 (θ̇) は 2 乗和多項式，
S4 (θ) は 2 乗和多項式(11)を考えることができる．実際，任意の (θ, θ̇) ∈ Θ に対して十分条件 (11) 式の行列はすべて半正定である．また，有界閉集合 (10) 式から r1 (θ) ≥ 0，r2 (θ̇) ≥ 0 であり，2 乗和(∀(θ, θ̇) ∈ Θ)． 多項式の性質 (9) 式から S1 (θ) ⪰ O，S2 (θ) ⪰ O，S3 (θ̇) ⪰(7) O，S4 (θ) ⪰ O であるため，r1 (θ)S1 (θ) ⪰ O，r1 (θ)S2 (θ)+
r2 (θ̇)S3 (θ̇) ⪰ O，r1 (θ)S4 (θ) ⪰ O が成り立つ．さらに，ここで，X(θ) は θ の 1 次関数 X0 + θX1 ，Ẋ(θ) は θ̇X1 ， M1 (ξ, θ) − r1 (θ)S1 (θ) ⪰ O，M2 (ξ, θ, θ̇) − r1 (θ)S2 (θ) −F (θ) は θ の 1 次関数 F0 + θF1 ，ξ は決定変数 ξ = r2 (θ̇)S3 (θ̇) ⪰ O，M3 (ξ, θ) − r1 (θ)S4 (θ) ⪰ O であるた
(X0 , X1 , F0 , F1 , γ)，Φ(θ) は Φ(θ) = AX(θ) + B2 (θ)F (θ) め，M1 (ξ, θ) ⪰ O，M2 (ξ, θ, θ̇) ⪰ O，M3 (ξ, θ) ⪰ O となであり，∗ は行列全体が対称になるように補完する．この る．したがって，十分条件 (11) 式をロバスト SDP(7) 式
とき，有解条件 (5) 式を達成するコントローラ (6) 式のゲ の PDLMI 制約条件と置き換えることにより，もとの問インは題の近似問題を作ることができる．この近似問題の実行可能解は必ずもとの問題の実行可能解になる．また，多K(θ) = F (θ)X(θ)−1項式行列 S1 (θ)，S2 (θ)，S3 (θ̇)，S4 (θ) の次数を高くする−1= (F0 + θF1 )(X0 + θX1 )(8) ことにより，十分条件 (11) 式は必要条件に近づき，近似の精度が高まる．しかし，2 乗和多項式に基づく手法は，によって与えられる．
θ の次元が高い場合や多項式行列 S1 (θ)，S2 (θ)，S3 (θ̇)，
ロバスト SDP(7) 式を解くには，PDLMI 制約条件がすS4 (θ) の次数が高い場合，最終的に解くべき問題の計算べての (θ, θ̇) ∈ Θ に対して満たされる必要があり，無限量が指数関数的に大きくなってしまう．
本の LMI を解かなければならないため，そのまま解くことは難しい．そこで，2 乗和多項式に基づきロバスト 3.5 他の方法による制御系の設計
SDP(7) 式を満たす近似解を求めることを考える．
2 乗和多項式に基づく制御の性能と比較するため，以下

3.32 乗和多項式多項式行列 S(θ) がある多項式行列 T (θ) を用いてS(θ) = T (θ)T T (θ)(9)と表すことができるとき，S(θ) を 2 乗和多項式という．このとき，S(θ) は対称行列，S(θ) の次数は偶数であり S(θ)
が 2d 次ならば T (θ) は d 次，任意の θ に対して S(θ) ⪰ Oであるという性質をもつ．の 2 つの方法による制御系の設計をする．
1 つ目の方法では，制御入力を力入力 uℓ = Cv (ẋw − ẋb )と仮定して状態空間表現を導出する．ここで，この状態空間表現は，アクティブサスペンション (サスペンション
外部の動力源から減衰力に相当する制御力を付加することによりバネ上 (車体) の振動を抑制する方法) を状態空間表現したものと等しい．この状態空間表現に対して，有解条件 (5) 式を達成するコントローラのゲイン (8) 式を求める．そして，実際の入力では Cv を Cv = ẋwu−ℓẋb によって決定する．この方法を以下，線形制御を適用した場合と呼ぶ．2 つ目の方法では，パラメータを θ とした線形パラメータ変動 (LPV) モデルを導出する．文献 [1] に従い，この LPV モデルに対して，ロバスト SDP(7) 式の実行可能解を求めることを考える．ロバスト SDP(7) 式の PDLMI
制約条件 M1 (ξ, θ)，M2 (ξ, θ, θ̇)，M3 (ξ, θ) はそれぞれ，θ，θ̇，θ に関して 1 次であり，さらに，M2 (ξ, θ, θ̇) は θ に関して 2 次であるため，θ2 を他とは独立なパラメータと見なし，θ，θ̇ の端点条件で置き換えることにより，有限個の LMI 制約条件へ帰着することができる．このように，
ロバスト SDP(7) 式の PDLMI 制約条件を有限個の LMI制約条件へ置き換えることにより，有解条件 (5) 式を達成するコントローラのゲイン (8) 式を求める．この方法を以下，従来法に基づく制御を適用した場合と呼ぶ．

4設計結果

表2

従来法に基づく制御を適用した場合

degXdegF

011101

表3γ

SeDuMitime[s]66.1106
66.1106
66.1106

0.644
0.640
0.682

γ

SDPT3time[s]66.1106
66.1106
66.1106

0.384
0.374
0.564

2 乗和多項式に基づく制御を適用した場合(a) X(θ) の次数が 0，F (θ) の次数が 1 のとき

degS1 ,degS2 ,degS4 degS32468————γ

SeDuMitime[s]γ

SDPT3time[s]66.0839 10.983 66.1106 7.786
66.1031 58.797 66.1106 40.749
66.0776 154.378 66.1106 48.25266.1056 348.459 66.1106 88.046

2 乗和多項式に基づきロバスト SDP(7) 式を満たす近似(b) X(θ) の次数が 1，F (θ) の次数が 0 のとき解を求める場合，解くべき問題の計算量が大きくなってし
SeDuMi
SDPT3
まうため，使用する計算機の性能により扱える問題の大degS1 ,degS2 ,degS4 degS3
γtime[s]γtime[s]きさに限界がある．本論文では，CPU として Intel Core22 12.6688 32.900 60.1324 21.372
i5 M560(2.67GHz)，メモリとして 3GB を搭載する計算
44 13.3362 159.833 58.8213 114.504
機上で，数値解析ソフトウェアとして MATLAB(VersionR2013a)，LMI パーサとして YALMIP(Version 19-Sep(c) X(θ) の次数が 1，F (θ) の次数が 1 のとき
2015)，SDP ソルバとして SeDuMi(Version 1.3.2) とSeDuMi
SDPT3
SDPT3(Version 4.0) を使用した．
degS1 ,degS2 ,degS4 degS3
γtime[s]γtime[s]各制御方法において，Q を単位行列，R を 1 とし，Se22 11.8799 30.561 59.5800 20.343
DuMi と SDPT3 それぞれで数値解を求めた．
44 13.3220 159.291 58.9789 112.712
線形制御を適用した場合では，SeDuMi と SDPT3 それぞれで γ = 66.1106 が得られた．従来法に基づく制御を 5 数値シミュレーション適用した場合の計算結果を表 2 に示す．ここで，X(θ) やF (θ) が定数行列 (X1 = O，F1 = O) であるとき X(θ) の次 5.1 数値シミュレーションの設定
一般路では，ストローク変位 0.005m 以下，ストローク数は 0，X(θ) や F (θ) が θ の 1 次関数 (X1 ̸= O，F1 ̸= O)
速度0.01m/s〜0.5m/s の領域が多用される．また，乗員であるとき X(θ) の次数は 1 である．表 2 より，X(θ) の次数や F (θ) の次数に関わらず，SeDuMi と SDPT3 それ は 0.5Hz から振動を感じ，特に，10Hz〜30Hz を不快にぞれで γ = 66.1106 が得られた．よって，これら 2 つの 感じる．そのため，自動車の乗り心地を評価するときこ制御方法では，SeDuMi と SDPT3 それぞれの数値解の れらの範囲が問題視される．よって，本論文では，これらの範囲の制振性能を保証することを考え，路面変位 xr
信頼性は同等であることが分かる．
2 乗和多項式に基づく制御を適用した場合の計算結果を にチャ―プ信号 (振幅 0.005m，時間に比例して周波数が表 3 に示す．ここで，多項式行列 S3 (θ̇) は θ̇ に関する 2 0.5Hz から 30Hz へ変化する正弦波) を加え，数値シミュ乗和多項式であり，X(θ) の次数が 0 であるとき多項式行 レーションをし，4 章の各制御方法を評価する．ここで，
列 S3 (θ̇) は存在しないため，「—」と表記する．また，本 各制御方法において，制御入力の大きさを調整するため，論文で使用した計算機の性能では，表 3 で扱った問題よ Q と R は試行錯誤により決定した．また，従来法に基づりも計算量の大きい問題の数値解を求めることはできな く制御を適用した場合と 2 乗和多項式に基づく制御を適かった．表 3 に見るように，SeDuMi と SDPT3 で異なる 用した場合の設計仕様を同等とするため，X(θ) の次数をγ の値が得られた．SeDuMi を使用して解いた計算結果で 0，F (θ) の次数を 1，多項式行列 S1 (θ)，S2 (θ)，S4 (θ) のは，多項式行列 S1 (θ)，S2 (θ)，S3 (θ̇)，S4 (θ) の次数を高 次数を 2，多項式行列 S3 (θ̇) の次数を 0 とした．
可変減衰係数 Cv を最小値 −500 に固定した場合 (ソフくしたとき，γ の値が増加する結果 (赤字部分) が得られた．これは 3.4 節で記した理論に反する結果である．一 ト固定) と最大値 500 に固定した場合 (ハード固定) にお方，SDPT3 を使用して解いた計算結果では，多項式行列 いて，縦軸を要求減衰力 N (= θ(Cs + Cv ))，横軸をストS1 (θ)，S2 (θ)，S3 (θ̇)，S4 (θ) の次数を高くしたとき，3.4 ローク速度 θ として，両者の関係を描いたものを図 2 に節で記した理論の通りに γ の値が減少する結果が得られ 示す．図 2 の赤線は減衰係数 Cv を調整することにより減た．よって，この制御方法では，SDPT3 の数値解の信頼 衰力可変式ショックアブソーバが減衰力を発生できる限性が高いことが分かる．5 章の数値シミュレーションでは， 界であり，コントローラ (6) 式から生じる要求減衰力はSDPT3 を使用して得たコントローラ (6) 式を使用する． 赤線で打ち切られて赤線の内側の領域で実現される．

5.2

数値シミュレーションの結果と考察

線形制御を適用した場合の状態の時間変化を，縦軸を要求減衰力 N ，横軸をストローク速度 θ として，描いたものを図 3 に示す．図 3 より，線形制御を適用した場合
では，コントローラ (6) 式から生じる要求減衰力が赤線の外側の領域にも存在していることが分かる．赤線の外側の領域にある要求減衰力は減衰力可変式ショックアブソーバがそのままの状態で実現することができないため，赤線で打ち切る．このとき，減衰力の不連続な切替が発生するため，乗員に制御の違和感を与えてしまう．
従来法に基づく制御を適用した場合と 2 乗和多項式に基づく制御を適用した場合の状態の時間変化を，縦軸を要求減衰力 N ，横軸をストローク速度 θ として，描いたものを図 4 と図 5 に示す．また，縦軸をバネ上 (車体) 変
位 xb ，ストローク速度 θ，可変減衰係数 Cv ，横軸を時間として，それらの時間変化を描いたものを図 6〜図 11 に示す．図 4 と図 5 より，これらの 2 つの制御方法では，コントローラ (6) 式から生じる要求減衰力が赤線の内側の領域に制限できていることが分かる．よって，減衰力可変式ショックアブソーバが要求減衰力をそのままの状態で実現することができるため，減衰力の不連続な切替が発生しない．図 10 と図 11 より，2 乗和多項式に基づく制御を適用した場合では，従来法に基づく制御を適用した場合と比べ，可変減衰係数 Cv を調整できる範囲を広く使うことができており，この意味で 2 乗和多項式に基づく手法は減衰力可変式ショックアブソーバの性能を十分に引き出しているといえる．

6

図2N vs θ の限界

図3線形制御：N vs θ

図4従来法：N vs θ

図52 乗和多項式：N vs θ

おわりに

本論文では，双線形システムとなるセミアクティブサスペンションを扱い，2 乗和多項式に基づく手法を適用した．4 章では，SDPT3 の数値解の信頼性が高いことを示した．また，他の方法に比べ γ の値が減少することを示した．5 章では，セミアクティブサスペンションへの有用性を示した．今後の課題は，状態や入力の制限を考慮した制御系の設計をすることである．

図6従来法：xb

図72 乗和多項式：xb

図8従来法：θ

図92 乗和多項式：θ

図 10従来法：Cv

図 112 乗和多項式：Cv

参考文献[1] R.Watanabe，K.Uchida and M.Fujita：A New LMI Approach to Analysis of Linear Systems with Scheduling
Parameter - Reduction to Finite Number of LMI Conditions-，Proceedings of the 35th IEEE CDC, 1663/1665(1996)[2] 市原裕之：二乗和に基づく制御系解析・設計，システム/制御/情報, Vol.55, No.5, 174/180 （2011）[3] 大石泰章：2 乗和多項式とその非線形制御への応用，システ
ム/制御/情報, Vol.58, No.11, 449/455 （2014）[4] S.Prajna，A.Papachristodoulou and F.Wu：Nonlinear
Control Synthesis by Sum of Squares Optimization：A
Lyapunov-based Approach，Proceedings of the 5th AsianControl Conference, 157/165 (2004)[5] Y.Fujisaki and Y.Oishi：Guaranteed cost regulator design：A probabilistic solution and a randomized algorithm，Automatica, Vol.43, No.2, 317/324 (2007)[6] 三平，大作，上村：非線形 H∞ 制御理論の限界と可能性―
セミアクティブサスペンションへの応用，システム/制御/情報, Vol.43, No.10, 544/552 （1999）