fc2ブログ

制御工学その29

8.1.3 安定性

状態空間表現から伝達関数表現への式は

P(s) = Y(s)/U(s) = C(sI-A)^-1*B+D

であった。このことから、Aの固有値は

|p*I - A| = 0

を満たすp_iであり、これは伝達関数P(s)の極に等しい。
ここから、

p_iの実部が全て負ならば、その時に限ってシステムは安定である。


8.2 状態方程式の解と遷移行列



状態方程式

x'(t) = A*x(t) + B*u(t)

において、u(t)=0とした零入力システム

x'(t) = A*x(t)

を考える。初期状態x(0)に対して、零入力システムの解は

x(t) = e^At*x(0)

である。ここでA*tは遷移行列と呼ばれ、

e^At = I + tA + t^2/2!*A^2 ... + t^k/k! * A^k ...

であらわされる無限級数である。

すると状態方程式の解は

x(t) = e^At*x(0) + int[0-t]{e^At*Bu(τ)}dτ

となる。

8.2.2 遷移行列の計算

今まででわかったように、システムの安定性は遷移行列e^Atに大きく依存する。
このe^Atを定義から求めるのは困難なので、ラプラス変換を利用して求めることが多い。

x'(t) = A*x(t)

をラプラス変換すると

s*X(s) - X(0) = A*X(s)

つまり

X(s) = (sI-A)^-1*x(0)

となり、これを逆ラプラス変換すると

x(t) = L^-1{(sI-A)^-1]*x(0)

となるので、

e^At = L^-1{(sI-A)^1}

となる。ここで、Aの固有値の実部が負であるなら、システムは任意の初期状態x(0)に対して
t→∞のときx→0となるため安定な振る舞いをする。特に、t→∞でx→0となるとき、漸近安定という。

例8.4 遷移行列の計算


A = |0,1|
|2,2|

(sI-A)^-1 = |s, -1|^-1
|2,s+2|

= 1/(s^2+2s+2) * |s+2,1|
| -2,s|

= 1/((s+1)^2+1^2) * {(s+1)*|1,0| + 1*| 1,1|}
|0,1} |-2,1|

ここで逆ラプラス変換して

e^At = e^-t*{|1,0|*cost + | 1,1|*sint}
|0,1| |-2,1}


となる。したがってx(t)は式

x'(t) = e^At*x(0)

つまり



|x1(t)| = e^-t*{|1,0|*cost + | 1,1|*sint} * |x1(0)|
|x2(t)| |0,1| |-2,1} |x2(0)|

= e^-t * | (cost+sint)*x1(0)+sint*x2(0)|
|-2sint*x1(0)+(cost-sint)*x2(0)|

となる。

問題8.2


A = | 0, 1|
|-6,-5|

のときの固有値e^Atを求める。


(sI-A)^-1 = (|s,0| - | 0, 1|)^-1 = |s, -1|^-1
|0,s| |-6.-5{ |6,s+5|

= 1/(s^2+5s+6)*|s+5,1| = 1/(s+2)(s+3) * |s+5,1| = k1/(s+2) + k2/(s+3)
| -6,s| | -6,s|



部分分数計算をする


k1(s+3) + k2(s+2) = |s+5,1|
| -6,s|

= (k1+k2)s + 3k1+2k2 = |s+5,1| = |s,0| + | 5,1|
| -6,s| |0,s| |-6,0|

k1+k2 = |1,0|
|0,1|

3k1+2k2 = | 5,1|
|-6,0|

3(I-K2) + 2k2 = | 5,1|
|-6,0|

-k2 = | 5,1| - 3I
|-6,0|

k2 = -(| 5,1| -3I) = -| 2, 1| = |-2,-1|
|-6,0| |-6.-3| | 6. 3|

k1 = I-k2 = | 3, 1|
|-6,-2|



よって、


e^-At = e^-2t*| 3. 1| + e^-3t*|-2,-1|
|-6,-2| | 6. 3|


となる。

8.3 制御系設計



8.3.1 可制御性と可観測性

制御の目的は、制御対象の状態x(t)を初期値x(0)から任意の状態に移すことである。
このような操作量は存在するだろうか。
任意の目標に到達できる操作量u(t)が存在することを、システムが可制御であるといい、システムが可制御であるかどうかは

可制御行列

V_c = {B,AB,...A^n-1*B}

が正則、rankV_cがnであることと等しい。

また、状態空間表現をベースに設計されるコントローラは通常状態変数x(0)をフィードバックする構造を持つ。
センサによってすべての成分が検出されるなら問題ないが、そうでないならば何らかの方法で状態変数を推測する必要がある。

制御量y(t)が検出可能である場合、y(t)は観測量と呼ばれるが、観測量と操作量から状態変数を正確に知ることが出来るならば
システムは可観測であるという。

システムが可観測であるかどうかは

可観測行列


V_o = | C|
| CA|
| .|
| .|
|C*A^n-1|


が正則である、すなわちrankV_o=nであることと等しい。


8.3.2 状態フィードバック

ここでは、状態変数x(t)が何らかの方法で利用可能であるとし、コントローラ

u(t) = K*x(t) + H*r

を用いることを考える。
ここで、K*x(t)は状態フィードバック、H*rは目標値rからのフィードバックである。
さらに、Kは状態フィードバックゲインと呼ばれる1*nのベクトル、Hはスカラーのフィードフォワードゲイン、rは定値の目標値である。

例8.5 垂直駆動アームのP-D制御

以前用いた垂直駆動アームのP-D制御コントローラは

u(t) = kp*e(t)-kD*y'(t)
e(t) = r - y(t)

であった。ただし、y(t) = θ(t),e(t) = r(t)である。
ここで、状態変数x1(t) = θ(t),x2(t) = θ'(t)とすると


u(t) = kp*(r-x1(t)) - kD*x2(t) = |-kp,-kD| * |x1(t)| + kp*r
|x2(t)|


となり、節頭で示したコントローラと同じ形式になる。

特に目標値がr=0であるとし、H=0とすると

u(t) = K*x(t)

となる。これは状態変数x(t)をフィードバックしているので状態フィードバックと呼ばれる。
状態フィードバックを用いたとき、状態方程式は

x'(t) = (A+BK)*x(t)

となるため、A+BKの固有値の実部がすべて負であるようにフィードバックゲインKが選ばれていると
システムは漸近安定となる。

まとめ

8.1.3

状態方程式

x'(t) = A*x(t) + B*u(t)

においてAの固有値

|p*I - A| = 0

を満たすp_iの実部が全て負ならば、その時に限ってシステムは安定である。


8.2

状態方程式の解は

x(t) = e^At*x(0) + int[0-t]{e^At*Bu(τ)}dτ

によってあらわされ、このときの

e^Atを遷移行列と呼ぶ。

遷移行列は主にラプラス変換を用いて計算することができる。

8.3

システムが実現したい目標値を達成するような操作量が存在するとき、システムは可制御であるといい、システムが可制御かどうかは
可制御行列が正則であるかどうかと等しい。

また、システムの制御量から内部状態を正確に観測することが出来るときシステムは可観測であるといい、システムが可観測かどうかは
可観測行列が正則であるかどうかと等しい。

8.3.2

状態フィードバックゲインと呼ばれるKを用いて

u(t) = K*x(t) + H*r

というコントローラを用いることが出来る。ここでr=0とし、H=0であるとすると
x'(t) = (A+BK)*x(t)となるので、システムが漸近安定かどうかは、

A+BKの固有値の実部がすべて負であるかを調べればよい。

スポンサーサイト



制御工学その28

8 さらに制御工学を学ぶ人のために




8.1 状態空間表現と安定性




8.1.1 状態空間表現

伝達関数表現は、システムの入出力関係のみに注目しているため、システムの内部信号がどのような振る舞いをしているのかを考慮できない。
さらに、初期の信号がすべて0であることを前提にしているため、初期値を考慮することが出来ない。
そこで、状態空間表現を用いる。

状態空間表現は、状態変数と呼ばれるシステムの内部信号x(t)を用いて、状態方程式と呼ばれる一階の微分方程式と出力方程式と呼ばれる
代数方程式とで記述される。
特に、システムが線形微分方程式で示される場合、状態空間表現は

状態方程式

x'(t) = A*x(t) + B*u(t)

出力方程式

y(t) = C*x(t) + D*u(t)

となる。
ただし、ここで

Aはn*nの行列
Bはn*1のベクトル
Cは1*nのベクトル
Dはスカラー

である。

このように、u(t),y(t)が共にスカラーであるようなシステムを1入出力系と呼ぶ。
また、システムが真にプロパーである場合、D=0である。

例8.1 垂直駆動アームの状態空間表現

いつもの垂直駆動アームは、運動方程式で表すと

Jθ''(t) = -cθ'(t) - Mlgθ(t) + τ(t)

となる。ここで、状態変数を


x(t) = |x1(t)| = | θ(t)|
|x2(t)| |θ'(t)|


u(t) = τ(t)

y(t) = θ(t)

と選ぶと、

x1'(t) = θ'(t) = x2(t)
x2'(t) = θ''(t) = 1/J * (-cθ'(t)-Mlgθ(t)+τ(t))
  = -Mlg/J*x1(t) -c/J*x2(t) + 1/J*u(t)

y(t) = θ(t) = x1(t)

となる。
そうすると、状態空間表現式を満たすには



x'(t) =|x1'(t)| = | x2(t)| = A * |x1(t)| + B*u(t)
|x2'(t)| |θ''(t)| |x2(t)|

y(t) = x1(t) = C*|x1(t)| + D*u(t)
|x2(t)|

を満たす必要があるため、

A = | 0, 1|
|-Mlg/J,-c/J|

B = | 0|
|1/J|

C = |1,0|

D = 0



となる。

8.1.2 伝達関数表現との関係

状態空間表現は伝達関数表現に容易に変換することができる。
初期状態をx(0) = 0とし、伝達関数表現の式をラプラス変換すると

sX(s) = AX(s) + BU(s)
よって

X(s) = (s*I-A)^-1*BU(s)

さらに

Y(s) = CX(s) + DU(s) = {C(sI-A)^-1*B+D}*U(s)

となるので、伝達関数P(s)は

P(s) = Y(s)/U(s) = C(sI-A)^-1*B+D

となる。


例8.2

さっき求めた垂直駆動アームの話

P(s) = Y(s)/U(s) = C(sI-A)^-1*B+D

に代入すると


P(s) = |1,0|*(|s,0| - | 0, 1|)^-1* | 0|
|0,s| |-Mlg/J,-c/J| |1/J|

= |1,0|*(| s ,-1|)^-1* | 0|
|Mlg/J,s+c/J| |1/J|

= |1,0|*1/(s^2+cs/J+Mlg/J)*(| s+c/J,1|)* | 0|
|-Mlg/J,s| |1/J|



= 1/(s^2+cs/J+Mlg/J)|s+c/J,1|*| 0| = 1/(s^2+cs/J+Mlg/J) * 1/J
|1/J|


となり、これは伝達関数表現と一致する。

一方、伝達関数表現から状態空間表現への変換は無数にあるため、よく知られたアルゴリズムがある。


例8.3 伝達関数表現から状態空間表現への変換

伝達関数表現

Y(s) = P(s)U(s)

P(s) = N(s)/D(s) = b2*s^2 + b1*s^2 + b0 / s^3+a2*s^2+a1*s+a0

というシステムが与えられたとき、状態空間表現への変換は、

V(s) = 1/D(s) * U(s)になるようなv(t)を考えれば

y(t) = b2*v''(t) + b1*v'(t) + b0*v(t)

u(t) = v'''(t) + a2*v''(t) + a1*v'(t) + a0*v(t)

と等価である。

したがって、状態変数を、

x1(t) = v(t)
x2(t) = v'(t)
x3(t) = v''(t)

と選ぶと、可制御標準形と呼ばれる、以下の状態空間表現が得られる。

可制御状態空間表現への変換


|x1'(t)| = |0,0,-a0|*|x1(t)| + |0|*u(t)
|x2'(t)| |1,0,-a1| |x2(t)| |0|
|x3'(t)| |0,1,-a2| |x3(t)| |1|

y(t) = |b0,b1,b2|*|x1(t)|
|x2(t)|
|x3(t)|



また、アルゴリズムは省略するが、可観測標準系と呼ばれる以下の状態空間表現に変換することも多い。

|x1'(t)| = |0,0,-a0|*|x1(t)| + |b0|*u(t)
|x2'(t)| |1,0,-a1| |x2(t)| |b1|
|x3'(t)| |0,1,-a2| |x3(t)| |b2|

y(t) = |0,0,1|*|x1(t)|
|x2(t)|
|x3(t)|



問題8.1

システムを状態空間表現で表す。
(1)伝達関数は

P(s) = 1/Ms^2+cs

正則表現化して

P(s) = 1/M / (s^2 + c/M*s)

a1 = c/M
a0 = 0
b0 = 1/M

x = x1=y(t),x2=y'(t)

とすると


|x1'(t)| = |0, 1|*|x1(t)| + |0|*u(t)
|x2'(t)| |0,-c/M| |x2(t)| |1|

y(t) = |1/M,0|*|x1(t)|
|x2(t)|


と言いたいけれど、どうやら違う形にして欲しいらしい。ぐぬぬ。


(2)

システムの微分方程式は

CL*y''(t)+RC*y'(t) + y(t) = u(t)

ここで
x1 = y(t)
x2 = y'(t)

とすると


|x1'(t)| = | x2(t)| = A * |x1(t)| + B*u(t)
|x2'(t)| = |-R/L*x2(t)-1/CL*x1(t)+1/CL*u(t)| |x2(t)|


y(t) = x1(t) = C*|x1(t)| + D*u(t)
|x2(t)|

A = | 0, 1|
|-1/CL,-R/L|

B = | 0|
|1/CL|

C = |1,0|

D = 0

またちょっと違うな。なんでだろう。

さらに伝達関数へ変換する。

公式は

P(s) = C(sI-A)^-1*B+D

なので



P(s) = |1,0| * (|s,0| - | 0, 1|)^-1*| 0|
|0,s| |-1/CL,-R/L| |1/CL|

(|s,0| - | 0, 1|)^-1 = | s, -1|^-1 = 1/(s^2+R/L*s+1/CL) * |s+R/L,1|
|0,s| |-1/CL,-R/L| |1/CL,s+R/L| |-1/CL,s|

P(s) = 1/(s^2+R/L*s+1/CL) * |1,0| * |s+R/L,1| * | 0|
|-1/CL,s| |1/CL|

P(s) = 1/(s^2+R/L*s+1/CL) * |s+R/L,1| * | 0| = 1/(CLs^2+CRs+1)
|1/CL|


となる。これはあってる。


まとめ

状態空間表現というものを考えた。
システムを状態変数と呼ばれる内部信号によって表現する方法なので、システムのそれぞれの内部状態が一手に分かる優れもの。
特に、システムが微分方程式で表されるならば、

状態方程式

x'(t) = A*x(t) + B*u(t)

出力方程式

y(t) = C*x(t) + D*u(t)

となる。
ただし、ここで

Aはn*nの行列
Bはn*1のベクトル
Cは1*nのベクトル
Dはスカラー


の2つ弐式で表すことが出来る。
また、この状態空間表現から伝達関数表現に変換することが出来る。
逆の場合、伝達関数から状態空間表現は無数の変換が考えられ(問題の時に示した)、簡単なアルゴリズムがよく知られている。


制御工学その27

7.3 フィードバック制御と周波数整形



7.3.1 フィードバック特性

ナイキストの安定判別により安定性を周波数領域で考えることが出来たが、
以下のフィードバック特性も周波数領域で考えることが出来る。

(a) 感度特性

制御対象の伝達関数P(s)がパラメータ変動や同定誤差の影響で実際には

P'(s) = P(s) + ΔP(s)

であったとすると、目標値r(s)から操作量y(s)への伝達関数

G_yr(s) = P(s)C(s)/1+P(s)C(s)

は実際には

G'_yr(s) = P'(s)C(s)/(1+P'(s)C(s)) = (P(s)+ΔP(s))C(s) / 1+(P(s)+ΔP(s))C(s)

となる。ここで、P(s)の変動に対するG_yrの感度をS(s)とし。

S(s) = (ΔG_yr(s)/G_yr(s))/(ΔP(s)/P(s))

ただし、ΔG_yr(s) = G'_yr(s) - G_yr(s)

と定義する。これはP(s)の変動率に対するG(s)の変動率の比である。
すると

S(s) = G'_yr(s)-G_yr(s)/G'_yr(s) * P'(s)/P'(s)-P(s) = (1- (P(s)C(s)/1+P(s)C(s))/(P'(s)C(s)/(1+P'(s)C(s)))) * P'(s)/P'(s)-P(s)

= (1 - P(s)*(1+P'(s)C(s))/(P'(s)*(1+P(s)C(s))) * P'(s)/(P'(s)-P(s)) = 1/1+P(s)C(s)

したがって、感度関数S(s)はP(s)が変動したときのG_yr(s)が受ける影響の指標になる。
目標値信号は通常、低周波成分を多く含むため、P(s)が変化したときのG_yr(s)が受ける影響を小さくするには、
|S(ωi)|を小さくする必要がある。


(b) 目標値追従特性

目標値r(s)から偏差e(s)への伝達関数は

G_er(s) = 1/1+P(s)C(s) = S(s)

である。目標値成分は低周波成分を多く含むため、目標値追従のためには|S(ωi)|を小さくする必要がある。


(c) 入力外乱抑制特性

入力外乱d(s)から制御量y(s)への伝達関数は

G_yd(s) = P(s)/1+P(s)C(s) = P(s)S(s)

である。入力外乱も低周波成分を多く含むため、入力外乱抑制のためには低周波領域で

|G_yd(ωi)| = |P(ωi)S(ωi)|を小さくすればいい。つまり、|S(ωi)|を小さくすれば、入力外乱を抑制することが出来る。


(d) 観測雑音除去特性

ノイズn(s)から制御量y(s)への伝達関数は

G_yn(s) = P(s)C(s)/1+P(s)C(s) = -T(s)

となる。ここで

T(s) = 1 - S(s)

であり、これを相補感度関数と呼ぶ。

ノイズは通常、高周波成分を多く持つため、観測雑音除去を行うためには|T(ωi)|を小さくする必要がある。



7.3.2 周波数整形

感度関数S(s)と相補感度関数T(s)には

S(s) + T(s) = 1

の関係が立つため、すべての周波数でS(s)とT(s)を同時に小さくすることは出来ない。
そこで、低周波領域でS(s)を小さく、高周波領域でT(s)を小さくするように、コントローラを設計すればいい。

また、低周波領域では |P(ωi)C(ωi)| >> 1であるため

S(ωi) = 1/1+P(ωi)C(ωi) ≒ 1/P(ωi)C(ωi)

と近似でき、高周波領域では |P(ωi)C(ωi)| << 1であるため、

T(ωi) = P(ωi)C(ωi)/1+P(ωi)C(ωi) ≒ P(ωi)C(ωi)

と近似することが出来る。

ここで、L(s) = P(s)C(s)の周波数整形の指標をまとめる。

(1)低周波領域

考慮する入力外乱と目標値の周波数がω_L以下であれば、ω_L以下で開ループ伝達関数のゲインを大きくし、感度特性、
目標値追従特性、入力外乱抑制特性を改善する。
また、ω→0におけるゲインの傾きが0なら定常偏差e_∞が残り、-20[dB/sec]で定常偏差e_∞が0になる。

(2)ゲイン交差周波数近辺

安定性を確保するため、ゲイン交差周波数ω_gc近辺のゲインの傾きを緩やかにし、十分な安定余裕をもたせる必要がある。
また、ω_gcを大きくすると速応性が向上し、ω_gcにおける位相余裕を大きくすると減衰性が向上する。

(3)高周波領域

考慮するノイズの周波数がω_H以上の高周波領域であるなら、ω_H以上で|L(ωi)|を小さくし、
観測雑音除去特性を改善する。


PIDコントローラを用いると、これらの指標をある程度実現できる。つまり、(1)を実現するにはC(s)に比例要素+積分要素を持たせればよく、
(2)を実現するにはC(s)に微分要素を持たせればよく、(3)を実現するにはC(s)に一次遅れなどのローパスフィルタを含ませれば良い。

まとめ

フィードバック特性を周波数領域で改善することを考えた。それぞれの入力は、それぞれ最も多く含む周波数領域をもつため、
それぞれが最も多い周波数領域でそれぞれ改善すればうまくいく。つまり

感度特性、目標値追従特性、入力外乱抑制特性はそれぞれの入力に低周波成分を多く含むため、低周波領域でこれらを改善できるようする。
観測雑音除去特性は、入力に高周波成分を多く含むため、高周波領域で改善できるようにする。

以下にそれぞれの指標をまとめる。

(1)低周波領域

考慮する入力外乱と目標値の周波数がω_L以下であれば、ω_L以下で開ループ伝達関数のゲインを大きくし、感度特性、
目標値追従特性、入力外乱抑制特性を改善する。
また、ω→0におけるゲインの傾きが0なら定常偏差e_∞が残り、-20[dB/sec]で定常偏差e_∞が0になる。

(2)ゲイン交差周波数近辺

安定性を確保するため、ゲイン交差周波数ω_gc近辺のゲインの傾きを緩やかにし、十分な安定余裕をもたせる必要がある。
また、ω_gcを大きくすると速応性が向上し、ω_gcにおける位相余裕を大きくすると減衰性が向上する。

(3)高周波領域

考慮するノイズの周波数がω_H以上の高周波領域であるなら、ω_H以上で|L(ωi)|を小さくし、
観測雑音除去特性を改善する。

そして、上に上げた3つは、PIDコントローラを用いればある程度実現することが出来る。
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR