目　　次 MATLAB/Simulink による制振入門講座

平成

18

3

22

-23

MATLAB/Simulink

福井大学工学部機械工学科　川谷　亮治

Email:[email protected] http://feedback.mech.fukui-u.ac.jp/

目　　次

1.

2.

2.1

2.2

2.3

2.4

3.

1

3.1

3.2

3.3

3.4

4.

4.1

4.2

4.3

4.4

4.5

4.6

5.

5.1

5.2

5.3

5.4

6.

6.1

6.2

H

6.3

6.4

H

6.5

LSDP

6.6

1

長大橋や超高層ビルなどの大型構造物は，その柔軟さゆえに地震や風などの外力を受け ることで振動が発生する。高速・軽量化された産業用ロボットアームや高所作業用ロング アームなどでも同様である。これらの対象において発生した振動を制御することは，居住 性・安全性・作業の高効率化などの点から必要な事項である。本稿では，フィードバック 制御を活用してこのような振動を制御する問題について考える。

本稿で具体例として取り上げる対象は，図

1.1

長大橋の簡易モデルとして考えることができる。

図

1.1:

第

2

設計において考慮すべきモデルのもつ不確かさについても触れる。本モデルは，厳密には 無限次元の振動モードの重ねあわせとして表すことができるが，そのうちの一つの振動 モードを取り上げ，それに対する解析とフィードバック制御による振動制御の考え方，設 計方法について第

3

4

(

)

5

2

MATLAB/SIMULINK

2

2.1

本稿における制御目的は 何らかの原因により柔軟ビーム中に発生した弾性振動をアク ティブに制御し，すばやく低減すること である。その目的を達成するために，柔軟ビー ムのある地点

(x = x

)

(x = x

)

f

(

1.1

)

2.2

今，時刻

t

x

y

(x, t)

EI ∂

y

∂x

+ m ∂

y

∂t

= 0 y

(0, t) = y

(L, t) = 0, ∂

y

∂x

¯ ¯

¯ ¯

= ∂

y

∂x

¯ ¯

¯ ¯

= 0

(2.1)

なお，式中の記号は以下のとおりである。

m :

L :

E :

I :

2

この偏微分方程式の解は

(

)

y

(x, t) = X

ϕ

(t) sin iπx

L (2.2)

上式は，柔軟ビームに発生した振動が，モード関数と呼ばれる

sin(iπx/L)

次に，外部から

(x

)

F(x, t)

EI ∂

y

∂x

+ m ∂

y

∂t

= F (x, t) (2.3)

上式に式

(2.2)

X

(m ϕ ¨

+ EI i

π

L

ϕ

) sin iπx

L = F (x, t)

その両辺に

sin(jπx/L)

[0, L]

sin

(m ϕ ¨

+ EI i

π

L

ϕ

) L

2 = Z

0

F (x, t) sin iπx

L dx (2.4)

を得る。本稿では，力

F(x, t)

x = x

f (t)

F (x, t) = f (t)δ(x − x

) (2.5)

と考えるので，これを式

(2.4)

i

sin(iπx/L)

ϕ

に関する微分方程式を得る。

mL

2 ϕ ¨

+ EIi

π

2L

ϕ

= sin iπx

L f (t) (i = 1, 2, · · · , ∞ ) (2.6)

ところで，実際の柔軟ビームでは，粘性摩擦等のエネルギ散逸が必ず存在するので，それ を考慮したのが次式である。

ˆ

a

ϕ ¨

+ ˆ b

ϕ ˙

+ ˆ c

ϕ

= ˆ d

f (i = 1, 2, · · · , ∞ ) (2.7)

ˆ

a

= mL

2 , ˆ b

= 2ζ

^p a ˆ

ˆ c

, ˆ c

= EIi

π

2L

, d ˆ

= sin iπx

L

2.3

2.3.1

現代制御理論では，制御対象に対して導かれた微分方程式を

1

(2.7)

x

x

=

⎡

⎢ ⎢

⎣ x

x

.. .

⎤

⎥ ⎥

⎦ ( x

=

"

ϕ

˙ ϕ

#

) (2.8)

x

i

u = f

diag { }

(

)

˙

x

= diag { A

, A

, · · ·} x

+

⎡

⎢ ⎢

⎣ b

b

.. .

⎤

⎥ ⎥

⎦ u (2.9)

ここで，

A

=

"

0 1

− ˆ c

/ˆ a

− ˆ b

/ˆ a

#

, b

=

"

0 d ˆ

/ˆ a

#

式

(2.9)

(2.8)

2.3.2

本稿で扱うのはフィードバック制御であるため，制御対象におけるどのような物理量が センサを利用して入手可能なのかを記述する必要がある。なぜなら，計測できないものを 利用して良好な特性をもつ制御が行えたとしてもそれは非現実的であるからである。

本稿で対象とする柔軟ビームの場合，図

1.1

x = x

h

(

)

y

= X

h 2 ∂

∂x

sin( iπx L )

¯ ¯

¯ ¯

x=xs

ϕ

(t)

= [c

c

· · · ]x

(2.10)

ここで，

c

= [ − h 2 ( iπ

L )

sin( iπx

L ) 0]

2.3.3

式

(2.9)

(2.10)

( x ˙

= A

x

+ b

u

y

= c

x

(2.11)

無限個の振動モードの重ね合わせとして柔軟ビームの振動に関する微分方程式

(2.7)

x

(2.11)

以下では，低次振動モードのみ

(

1 ∼ 3

)

( x ˙ = Ax + bu

y = cx (2.12)

ここで，

x =

⎡

⎢ ⎣ x

x

x

⎤

⎥ ⎦ , A =

⎡

⎢ ⎣

A

0 0 0 A

0 0 0 A

⎤

⎥ ⎦ , b =

⎡

⎢ ⎣ b

b

b

⎤

⎥ ⎦ , c = ^h c

c

c

i

この場合，無視した

4

(2.12)

5

2.4

それでは，状態空間モデル

(2.12)

表

2.1:

L 1.8 [m]

m 0.2421 [kg/m]

h 0.003 [m]

EI 5.063 [N m

] E × I

ζ

0.05

(1

)

ζ

0.008

(2

)

ζ

0.006

(3

)

x

3L/4 [m]

x

L/4 [m]

0.078 [N/V ]

1/(20 × 10

) [1/V ]

アクチュエータとセンサをそれぞれ両端から

L/4

4

(non-colocation)

(colocation)

例題

2.1

2.1

MATLAB

(3

)

bmodel.m

function bm=bmodel()

%

(3

)

%

L=1.8; m=0.2421; h=0.003; EI=5.063;

zeta1=0.05; zeta2=0.008; zeta3=0.006;

xa=3*L/4; xs=L/4;

ampgain=0.078; sensorgain=1/(20*1e-6);

a=m*L/2;

c1=EI*pi^4/2/L^3; d1=sin(pi*xa/L); b1=2*zeta1*sqrt(a*c1);

c2=EI*(2*pi)^4/2/L^3; d2=sin(2*pi*xa/L); b2=2*zeta2*sqrt(a*c2);

c3=EI*(3*pi)^4/2/L^3; d3=sin(3*pi*xa/L); b3=2*zeta3*sqrt(a*c3);

cc1=-h/2*(pi/L)^2*sin(pi*xs/L);

cc2=-h/2*(2*pi/L)^2*sin(2*pi*xs/L);

cc3=-h/2*(3*pi/L)^2*sin(3*pi*xs/L);

A1=[0,1;-c1/a -b1/a]; B1=[0;d1/a]; C1=[cc1 0];

A2=[0,1;-c2/a -b2/a]; B2=[0;d2/a]; C2=[cc2 0];

A3=[0,1;-c3/a -b3/a]; B3=[0;d3/a]; C3=[cc3 0];

z2=zeros(2);

A=[A1 z2 z2;z2 A2 z2;z2 z2 A3]; B=[B1;B2;B3]*ampgain;

C=[C1 C2 C3]*sensorgain; D=0;

bm=ss(A,B,C,D);

%EOF

例題

2.2

(impulse)

>> bm=bmodel;

>> t=0:0.01:10;

>> impulse(bm,t)

得られた図から，本対象は振動性が強く，しかもその減衰がゆるやかで，インパルス状の 外乱が加えられてから応答が収束するまでに

8[s] ∼ 10[s]

図

2.2:

3 1

3.1

フィードバック制御による振動制御の本質を理解する目的で，本章では次図に示す

1

k

μ

図

3.3: 1

本系に対する微分方程式は力の釣り合いから次式で与えられる。これは柔軟ビームにおけ る一つの振動モードを取り出した微分方程式と等価である。

M z ¨ + μ z ˙ + kz = u (3.13)

上式に対して，状態量を

x = [z, z] ˙

z

M = 1[kg]

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

˙ x =

"

0 1

− k − μ

# x +

"

0 1

# u y = ^h 1 0 ⁱ x

(3.14)

例題

3.1

k, μ(mu)

(3.14)

mck.m

function s1=mck(k,mu)

%

%

s1=ss([0,1;-k,-mu],[0;1],[1,0],0);

%EOF

3.2

演習

3.1

mck.m

k, μ

BODE

(bode)

>> t=0:0.01:10; %

>> w=logspace(-1,2,100); %

>> s1=mck(10,1); %

>> impulse(s1,t) %

>> bode(s1,w) % BODE

この演習を通して，以下のことがわかる。

1. k < 0

μ < 0

2. k > 0

μ = 0

3. k > 0

μ > 0

y = 0

1

(

)

3

(

)

4. μ

k

BODE

5.

k

μ

BODE

6. μ

ところで，これらの結果をより客観的に表現するために，次に示す代数方程式とその根

λ

(3.14)

A

(

)

A

s

+ μs + k = 0 (3.15)

これに対して以下に示す結果が得られる。

1.

(

)

2.

約

8/μ[s]

3.

以上より，振動特性を考える上で，特に

μ

3.3

それでは，式

(3.13)

(

)

k, μ

u = − k

z − k

z ˙ = − ^h k

k

ⁱ

"

z

˙ z

#

= − Kx (3.16)

というように操作量

u

K = [k

, k

]

上式を式

(3.14)

⎧ ⎪

⎪ ⎪

⎨

⎪ ⎪

⎪ ⎩

˙

x = (A − BK)x =

"

0 1

− (k + k

) − (μ + k

)

# x y = ^h 1 0 ⁱ x

(3.17)

となり，特性方程式は

s

+ (μ + k

)s + (k + k

) = 0 (3.18)

で与えられる。

k

, k

(3.17)

(

(3.14))

(

(3.17))

3.4

式

(3.14)

μ = 2, k = 101

− 1 ± 10j

(

3.4

)

1[s]

2

図

3.4:

>> t=0:0.01:10; %

>> s1=mck(101,2); %

>> eig(s1) %

>> impulse(s1,t) %

Case 1

まず，応答から振動そのものをなくしてしまうことを考えよう。

3.2

μ = 2

μ = 8

− 4 ± 4j

(s + 4 + 4j)(s + 4 − 4j) = s

+ 8s + 32 (3.19)

u = 69z − 6 z ˙ (3.20)

が望ましい制御則となる。

Case 2

次の考え方は，振動をなくすのではなく，その減衰特性を高めてやろうとするものであ る。整定時間は，

μ

μ

2

8

u = − 6 z ˙ (3.21)

で与えられる。

実際に数値計算により単位ステップ応答を求め，両者の比較を行う。

>> t=0:0.01:5; %

>> yopen=impulse(s1,t); %

>> sc1=mck(32,8); % Case 1

>> sc2=mck(101,8); % Case 2

>> ycl1=impulse(sc1,t); %

>> ycl2=impulse(sc2,t); %

>> plot(t,[ycl1,ycl2,yopen]), grid on %

図

3.5: 1

破線が

Case1

Case2

Case1

Case2

(

1[s])

ところで，フィードバック制御はシステムに対して操作量を与えることで実現するので，

Case1

Case2

>> ucl1=impulse([0,1;-32,-8],[0;1],[-69,6],0,1,t);

>> ucl2=impulse([0,1;-101,-8],[0;1],[0,6],0,1,t);

>> plot(t,[ucl1,ucl2]), grid on

図

3.6:

破線が

Case1

Case2

必要としていることがわかる。制御対象によっては，

Case1

減衰を高めることを目指すほうが現実的である。特に，大型構造物の場合，いかに少ない 操作量で希望する特性を実現するかは非常に重要なことである。

Case2

− 1 ± 10j

− 4 ± 9.22j

(

)

>> eig(sc2) % Case2

4

前章では，式

(3.14)

(3.16)

yes

4.1

最初に，可制御性 を定義する。対象とするのは次に示す状態方程式である。

˙

x = Ax + Bu (4.22)

可制御性 定義

4.1

x(0)

t

x(t

)

0

u(t) (0 ≤ t ≤ t

)

(4.22)

図

4.7:

なお，システムが可制御であることを対

(A, B)

(

)

ところで，定義自体では，与えられたシステムが可制御であるかどうかを判定するには 不便なので，その判定に関連した重要な定理を紹介する

(

)

可制御性の判定法

定理

4.1

(A, B)

V = [B AB A

B · · · A

B] (4.23)

(

)

rank(V ) = n

つまり，状態方程式が与えられたとき，可制御行列

V

例題

4.1

(2.12)

>> bm=bmodel; %

>> [A,B,C,D]=ssdata(bm); % A,B,C,D

>> V=ctrb(A,B); %

>> rank(V) %

4.2

可制御性と状態フィードバック制御

定理

4.2

(A, B)

A − BK

指定できる行列

K

状態方程式

(4.22)

u = − Kx

˙

x = (A − BK)x (4.24)

で与えられる。上の定理は状態フィードバック制御により閉ループ極を任意に指定可能で あることを意味している。

ところで，閉ループ極を任意に指定できる状態フィードバックゲイン

K

これが 極配置法 と呼ばれる設計法である。この設計法では，閉ループ極をどのように定 めるかが鍵となる。一般的な考慮すべき点を以下に列挙する。

1.

2.

3.

t

t

= 4

支配極^(∗)の実部の絶対値

4.

極の実部の絶対値

≥

5.

(*)

4.3

それでは，柔軟ビームに対して極配置法による状態フィードバックゲインの設計ならび にその評価を行う。ここでは，設計仕様として，単位インパルス外乱入力に対して，約

1[s]

程度で応答を整定させることを目指す。

演習

4.1

>> help place %

place

>> bm=bmodel; %

>> eig(bm) %

>> [A,B,C,D]=ssdata(bm); % A,B,C,D

>> t=0:0.01:5; %

>> w=logspace(-1,2,100); %

>> K=place(A,B,[???,???,???,???,???,???]); %

>> impulse(A-B*K,B,[C;-K],[D;0],1,t); %

>> bode(A-B*K,B,C,D,1,w); %

4.4

極はシステムの応答に強く関係している。極配置法は閉ループ極を直接指定できるので，

状態フィードバックゲインを設計するための有力な手法の一つとして挙げることができる。

ところで，前章の

1

(

)

4.4.1

˙

x = Ax + Bu (4.25)

に対して，次に示す 評価関数

J

u

J = Z

0

(x

Qx + u

Ru) dt (4.26)

ここで，行列

Q, R

(

)

(

)

Q

R

n

x

Q = diag { q

, q

, · · · , q

}

(4.26)

1

Z

0

x

Qx dt = Z

0

X

i=1

q

x

dt = X

i=1

q

Z

0

x

dt

となる。これより，評価関数の第

1

2

2

2

(4.26)

4.5

最適レギュレータ問題の解を求めるために，次式に示す

Riccati

A

X + XA − XBR

B

X + Q = 0 (4.27)

Riccati

X

dt d (x(t)

Xx(t)) = x ˙

Xx + x

X x ˙

= (Ax + Bu)

Xx + x

X(Ax + Bu)

= x

(A

X + XA)x + u

B

Xx + x

XBu

= − (x

Qx + u

Ru) + (u + R

B

Xx)

R(u + R

B

Xx) 2

(4.25)

4

Riccati

方程式

(4.27)

[0, ∞ )

J = − x

(t)Xx(t)

¯ ¯

¯

0

+ Z

0

(u + R

B

Xx)

R(u + R

B

Xx) dt

最適レギュレータ問題の解

u = − R

B

Xx (4.28)

と選ぶことが評価関数

J

2

u

J

J

= x

(0)Xx(0) (4.29)

で与えられる。なお，最適レギュレータ問題の設定は状態フィードバック制御を前提とし ていないが，結果として得られる最適操作量

u

K

状態フィードバックゲイン

K = R

B

X (4.30)

である状態フィードバック制御

(u = − Kx)

例題

4.2

3.4

1

Q, R

3.4

(

1[s])

Q =

"

1 0 0 60

#

, R = 1

− 4 ± 9.22j

>> help lqr %

lqr

>> s1=mck(101,2);

>> [A,B,C,D]=ssdata(s1);

>> t=0:0.01:5;

>> K=lqr(A,B,diag([1,60]),1); %

>> impulse(A-B*K,B,[C;-K],[D;0],1,t);

4.6

演習

4.2

1[s]

>> bm=bmodel;

>> [A,B,C,D]=ssdata(bm);

>> t=0:0.01:5;

>> Q=...; R=...; %

Q, R

>> K=lqr(A,B,Q,R); %

>> impulse(A-B*K,B,[C;-K],[D;0],1,t);

演習

4.3

4.2

1 ∼ 3

1

5

状態方程式はシステムの動特性を表現している。

˙

x = Ax + Bu (5.31)

これに対して状態フィードバック制御

u = − Kx

K

x

6

1

5.1

˙

x = Ax + Bu (5.32)

において，行列

A, B

u

ˆ ˙

x = Aˆ x + Bu (5.33)

は計算機等を利用して解くことができ，式

(5.33)

x ˆ

x ˆ

x

(

)

(5.32)

(5.33)

e = x − x ˆ

˙

e = Ae (5.34)

式

(5.32)

x

x(0)

し，微分方程式

(5.33)

x(0) ˆ

x(0) − x(0) = ˆ e(0) 6 = 0

x(0)

e(0)

e( ∞ ) = 0

lim

t→∞

x ˆ = x

(5.34)

A

lim

t→∞

x ˆ = x

x ˆ

x

no

しかし，式

(5.34)

e

e

ˆ

x

そこで，式

(5.33)

x ˆ

ˆ y = C x ˆ

を考える。これは，観測量

y

y − y ˆ

(5.33)

H

ˆ ˙

x = Aˆ x + Bu + H(y − C x) ˆ (5.35)

上式と式

(5.32)

˙

e = (A − HC)e (5.36)

となる。状態フィードバック制御

u = − Kx

A

A − BK

が，式

(5.35)

H

ることが，図

5.8

(5.34)

A

A − HC

図

5.8:

(5.35)

H

したがって，行列

H

A − HC

e

(5.35)

x ˆ

x

(5.35)

(

)

H

状態量が推定できるならば，それを利用して状態フィードバック制御を実現できる。つ まり，全状態オブザーバを用いた フィードバック制御器 は次式で与えられる。

( x ˆ ˙ = A x ˆ + Bu + H(y − C x) ˆ

u = − K x ˆ (5.37)

あるいは，第

2

u

1

u

y

(5.38)

y

u

( x ˆ ˙ = (A − BK − HC)ˆ x + Hy

u = − K x ˆ (5.38)

5.2

それでは，可観測性を定義しよう。

可観測性

定義

5.1 ₍

˙

x = Ax

y = Cx (5.39)

において，ある有限時刻

t

y(t) (0 ≤ t ≤ t

)

x(0)

(5.39)

なお，システムの可観測性は行列

A

C

(C, A)

可制御性のときと同様に，与えられたシステムが可観測であるかどうかを判定するため の重要な定理を紹介する。

可観測性の判定法

定理

5.1

(C, A)

W =

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎣ C CA CA

.. . CA

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎦

(5.40)

がフル

(

)

rank(W ) = n

前章で述べた可制御性は，任意の初期状態

x(0)

x = 0

u

x

x(0)

x

(C, A)

(A

, C

)

(A, B)

(B

, A

)

もう一つの重要な定理を紹介する。

可観測性とオブザーバゲイン

定理

5.2

(C, A)

A − HC

指定できる行列

H

　 つまり，

システムが可観測であるならば，適切に選定したオブザーバゲイン

H

(5.35)

x ˆ

x

ところで，オブザーバゲイン

H

(C, A)

(A

, C

)

A

− C

H ˆ

H ˆ

A − H ˆ

C

H = ˆ H

図

5.9: A − BK

A − HC

それでは，具体的にどこに配置するかということになるが，推定した状態量

x ˆ

(

A − BK

)

(

5.9

)

A − HC

例題

5.1

1

>> s1=mck(101,2); %

>> W=obsv(s1); %

>> rank(W) %

5.3

演習

5.1

>> bm=bmodel;

>> [A,B,C,D]=ssdata(bm);

>> t=0:0.01:5;

>> K=place(A,B,[??,??,??,??,??,??]); %

>> H=place(A’,C’,[??,??,??,??,??,??])’; %

>> ak=A-B*K-HC; bk=H; ck=-K; dk=0; %

>> [ac,bc,cc,dc]=feedback(A,B,C,D,ak,bk,ck,dk,1); %

>> ycl=impulse(ac,bc,cc,dc,1,t); %

>> plot(t,ycl), grid on %

5.4

可観測なシステムに対して，全状態オブザーバを利用することで状態量

x

x ˆ

u = − K x ˆ

e = x − x ˆ

u = − Kx

状態空間モデル

( x ˙ = Ax + Bu

y = Cx (5.41)

に対して制御器

(5.38)

( x ˙ = Ax − BK x ˆ ˆ ˙

x = HCx + (A − BK − HC)ˆ x (5.42)

これに対して状態量を

[ x

, (x − x) ˆ

]

"

˙ x

˙ x − x ˆ ˙

#

=

"

A − BK BK

0 A − HC

# "

x x − x ˆ

#

(5.43)

を得る。これより，閉ループシステムに対する状態方程式の極は，レギュレータ極

(A − BK

)

(A − HC

)

K

H

ˆ

x

u = − K x ˆ

>> eig(ac) %

6

本章では，

H

(Loop Shaping Design Pro- cedure:

LSDP )

6.1

4

5

3

2

ように，柔軟ビームに対する状態空間モデルは，発生する振動が無限個の振動モードの重 ね合わせとして表現されることに対応して，無限次元をもつ。したがって，設計した制御 器を実制御対象に適用した場合，無視した高次振動モードの影響で応答特性が悪化したり，

場合によっては不安定になってしまう

(

)

H

以降の節で

H

1

6.2 H

H

(

)

H

H

H

定義

6.1

P (s)

H

k P k

= sup

ω

| P(jω) | (6.44)

つまり，伝達関数

P (s)

BODE

H

例題

6.1

1

BODE

H

(norm)

>> s1=mck(101,2); %

>> [A,B,C,D]=ssdata(s1);

>> [num,den]=ss2tf(A,B,C,D,1); %

>> w=logspace(0,2,200);

>> [g,p]=bode(num,den,w); %

>> semilogx(w,20*log10(g)), grid on %

>> xlabel(’Freq[rad/s]’), ylabel(’Gain[dB]’)

>> max(g) %

>> norm(s1,inf) % Hinf

図

6.10: H

ω

k P k

≥ | P (jω) | (6.45)

6.3

ロバスト安定性に関連する重要な定理が スモールゲイン定理 である。図

6.11

図

6.11:

図中，

P (s)

∆(s)