第2章 数学モデル

第 2 章 数学モデル
2.1
システムの一般的性質
全く異なるシステムが, 同じ数式で表されることがしばしばある.
最初の例は, 抵抗である. 電気抵抗 R は, オームの法則で式 (2.1) のように表される
V = RI
(2.1)
ここで V は電圧 (抵抗を横切る駆動ポテンシャル), I は抵抗を流れる電流である. V が, 「横切る」
変数であり, 「作用力」の尺度であることに注意. 一方, I は, 「通す」変数であり, そして「流れ」
の尺度である.
一般的に, 「作用力」の変数を ψ, 流れの変数を ζ と置くと, オーム法則は, 式 (2.2) と書ける.
ψ = Rζ
(2.2)
ここで, R は一般化された抵抗である.
図 2.1 に異なるシステムの一般化された抵抗の概念の応用を示す.
• 機械的なダッシュポットでは, 力 F がプランジャーに加えられると, それは速さ v で動く. v
は F に比例する. F は ψ に, v は ζ に対応する. 比例定数 Rm は, ダッシュポット内の流体の
粘性に関係する. 機械的抵抗である. 実際に, Rm は, 衝撃吸収としてのダッシュポットの性
能を決め, もっと一般的には, 制動係数として知られている.
F = vRm
(2.3)
• 液体の流れでは, 一般化されたオームの法則は, Poiseuille の法則を仮定する. 硬い管を流れ
る流体の体積 Q は, 管の両端の圧力差 ∆P に比例する (管の半径の 4 乗に比例, 管の長さに反
比例). 流体の抵抗 Rf は, 流体の粘性と管の長さに比例し, 管の断面積の 2 乗に反比例する.
∆P = QRf
(2.4)
• 熱抵抗 Rt は, 材料の熱伝導率の逆数に比例する.
∆Q = QRt
(2.5)
• 化学システムでは, 濃度が違う 2 つの液体が透過性のある膜でしきられているとき, 膜を透過
するフラックス Q は, 濃度差 ∆φ に比例する. これはフィックの拡散の法則.
∆φ = QRc
ここで Rc は, 拡散抵抗.
10
(2.6)
図 2.1: 抵抗 (a) 力学系, (b) 流体, (c) 熱, (d) 化学系
11
2 つめの例は, 貯蔵 (storage) である.
• 電気系では, 単位電圧 (V ) あたりに蓄えられている電荷 (q) の量として定義されるキャパシ
タンスである.
q
(2.7)
C=
V
q は, コンデンサに流れ込む電流のよる全ての電荷の積分である.
t
Idt
(2.8)
q=
0
式 (2.7) および式 (2.8) から式 (2.9) を得る.
ψ=
1
C
t
ζdt
(2.9)
0
• 機械系では, コンプライアンス (弾性バネ) である. コンプライアンス Cm は, バネの伸び・縮
みの大きさ x を決める. この特性は, バネの硬さ・弾性率に反比例する. バネが軟らかい程,
与えられた力でよく伸びる.
x
=
x
=
F
=
F Cm
t
vdt
0
t
1
vdt
Cm 0
(2.10)
(2.11)
(2.12)
• 流体系では, 流体をみたした風船. コンプライアンス (単位圧力変化あたりの体積変化) は, 圧
力の変化に対する膨張・収縮することによる体積を決める. コンプライアンスは, 風船の材料
の弾性によって決まる. 風船の中の流体の圧縮性によるコンプライアンスへの寄与はとても
小さい.
∆V = ∆pCm
(2.13)
• 熱の系では, 温度差によってある媒質の中に蓄えられる熱の量として表される熱溜である (単
位温度差あたりの熱). 熱容量は大きさと媒質の熱に依存する.
q = ∆θCt
(2.14)
• 化学系では, 貯蔵は, 化学種が存在する流体の全体積によって表される. 例えば, 与えられた
体積にとって, 化学種の全質量は, 濃度に比例する. これは濃度の定義.
m = cv
(2.15)
ここで m は質量, c は濃度, v は体積.
抵抗はエネルギーを散逸, 貯蔵はエネルギーを蓄える. 電気では, 電圧 (V) と電流 (A) の積はパ
ワー (W). 機械では, 力 (N) と速度 (m/s) の積はパワー. 流体では, 圧力 (Pa) と流速 (体積流量率
m3 /s) の積がパワー.
どんなシステムでも, 「作用力」が抵抗を通して流れを作るために適用されると, エネルギーが
散逸する. 一方, 貯蔵要素では, 静的あるいはポテンシャルエネルギーが蓄積される. 機械系のばね
では, x の伸長を起こす力 F は, ポテンシャルエネルギー F x の保存となる.
12
図 2.2: Strage (a) 力学系, (b) 流体, (c) 熱
13
最後の性質は, 慣性である. 電気系における動的なエネルギーの貯蔵はインダクタンス L である.
インダクタンス L は, 与えられた電流の時間変化を生成するために必要な電圧として定義される.
V =L
dI
dt
(2.16)
一般化して,
dζ
(2.17)
dt
機械系では, ψ は F , ζ は速度で, 式 (2.17) はニュートンの第 2 法則. この場合には, 慣性は質量. 流
体では, 流体の加速度が圧力差に比例する. 熱や化学システムには, 動的なエネルギーの貯蔵は存
ψ=L
在しない.
2.2
システム要素の組合せモデル
基本的なモデル要素は, 線形であると仮定する. 実際には違う. 電気抵抗は電流が流れることで
発熱し, 電気抵抗が増加する. 流体では, 層流では, 流体抵抗は, 比較的に一定である. しかし, 流れ
が増加して乱流になると, 抵抗は流れそれ自身の関数になる. 非線形な作用力が増加すると, これ
らの異なるシステムの類似性が減少する. しかし, 簡単なモデルからスタートすることが賢い. 第
一近似として, 線形解析は, 生体システムのモデリングで重要である.
図 2.3: 一般化したシステム構成要素からなる単純なモデル
図 2.3 は, 3 つの基本要素が直列・並列に接続されている簡単な回路である. 「0」は, 「作用力」
ψ の基準である. 電気回路では, 接地と呼ばれる. ψa と ψb は, ノード a と b で, 「横切る変数」の
(0 に対する) 値である. ζ1 , ζ2 , ζ3 は, 要素 E1, E2, E3 を流れる, 「流れの変数」の値である. これ
らには, 2 つの基本的な式がある.
1. 任意の閉回路で, 「横切る変数」の代数和は零である.
(ψa − ψb ) + (ψb − 0) + (0 − ψa ) = 0
14
(2.18)
2. あるノードにはいってくる「流れの変数」の代数和は零である.
ζ1 + (−ζ2 ) + (−ζ3 ) = 0
(2.19)
電気回路では, キルヒホッフの法則. E1, E2, E3 は, 基本要素であるが, それ自身がネットワークで
あってもよい.
図 2.4: 抵抗と strage の要素の組合せ
図 2.4 上半分は, 抵抗と貯蔵の要素 (電気, 流体, 熱, 化学) の組合せのシステムの特性である.
抵抗とコンプライアンスが結合したものは 図 2.4 下半分に示すように, 少し違う.
(a) の場合について考える. 2 つのダッシュポットが並列のとき, 同じ速度 v で動かすための力 F
15
は,
F = v(Rm1 + Rm2 )
(2.20)
R = Rm1 + Rm2
(2.21)
になる. 組み合わせた抵抗 R は,
になる. 機械的な抵抗 2 つの並列は, 他のタイプの抵抗とは異なる.
(c) の場合に付いて考える. 2 つのコンプライアンスが Cm1 と Cm2 のバネが並列のとき, 同じ量
x の変位について,
x = F1 Cm1 = F2 Cm2
(2.22)
となる, F1 と F2 が全体の力なので, バネ全体のコンプライアンスは, C = x/F であるから,
1
1
1
=
+
C
Cm1
Cm2
(2.23)
となる. 他の場合のキャパシタンスの並列とは異なる.
機械系と電気系のモデルを変換するときには注意すること.
2.3
生体システムの線形モデル：2 例
図 2.5: 肺のモデル
1 つめは, 肺のモデル (図 2.5). 気道は 2 つに分類される. 大きい中央の気道と, 小さい周辺の気
道である. 流体の抵抗を Rc と Rp とおく. 肺胞に入った空気は, 胸腔の拡大 (体積ははいってきた空
気と同じ) を発生する. 直列につながった, 肺のコンプライアンス (CL ) と, 胸壁のコンプライアン
ス (CW ). 呼吸器にはいってきた空気の内, 一部は, 中央の気道のコンプライアンスとガスの圧縮性
により, 肺胞に入らない. 通常の呼吸周波数で, 正常な状態では, このわきにそれる分は少ない. し
かし, 病気が末梢の気道の閉塞 (Rp の増加) や, 肺や胸壁の硬化 (CL や CW の減少) を引き起こし
たら, 重要になる. このコンプライアンスを CS として, CL および CW に並列にする. 圧力を次の
ように置く. Pa0 は気道の入口, Paw は中央の気道, PA は肺胞, Ppl は胸膜 (肺の組織と胸壁の間),
P0 は周囲の圧力 (零と設定できる).
16
体積流量 Q の空気が入ると仮定する. Q と Pa0 の関係を求める.
ノード Paw に, キルヒホッフの第 2 法則を適用する. 肺胞にはいる空気を QA とする. わきにそ
れるのは, Q − QA になる.
1
∆P = QRf , ∆P =
Cm
Qdt
(2.24)
であることに注意して, CS , RP , CL と CW の閉回路にキルヒホッフの第一法則を適用する.
1
1
1
+
(2.25)
(Q − QA )dt
RP QA +
QA dt =
CL
CW
Cs
RC と CS を含む回路にキルヒホッフの第一法則を適用する.
1
(Q − QA )dt
Pa0 = RC Q +
CS
(2.26)
式 (2.25) および式 (2.26) を t に関して微分して, QA を消去して一つにまとめる.
QA
RP
+
dt
1
1
+
CL
CW
QA
dPa0
dt
dPa0
CS
dt
RP CS
=
=
QA
=
dQA
dt
=
1
(Q − QA )
CS
1
dQ
RC
+
(Q − QA )
dt
CS
dQ
CS RC
+ Q − QA
dt
dPa0
dQ
+ CS RC
+Q
−CS
dt
dt
d2 Pa0
d2 Q dQ
−CS
+
C
R
+
S
C
dt2
dt2
dt
(2.27)
(2.28)
(2.29)
(2.30)
(2.31)
dPa0
dQ
+ CS RC
+ Q (2.32)
−CS
dt
dt
1
dPa0
dQ
− CS RC
=
CS
(2.33)
CS
dt
dt
1
1 dPao
RC
1
1
dQ
1
d2 Q
= RC 2 +
+
+
+
+
Q
(2.34)
RP CT dt
dt
CS
RP CT
dt
RP CS CL
CW
d2 Pa0
d2 Q
dQ
+
+ CS RC RP 2 + RP
2
dt
dt
dt
d2 Pao
dt2
=
1
1
+
CL
CW
ここで, CT は, 次のように定義される.
CT =
1
1
1
+
+
CL
CW
CS
−1
(2.35)
2 つめの例は骨格筋 (図 2.6). F0 は収縮要素によって発生する力. F が実際の力. R は組織の粘
性. CP と CS は, 並列弾性要素と直列弾性要素. それぞれ, サルコメアと腱の貯蔵する性質を反映
している. もし, CP が x 伸ばされたら, R と CS 全体も同じだけ伸ばされる. 2 つの並列の要素の
力の和が F . CS の伸びを x1 とすると, R と F0 の伸びは x − x1 になる. ダッシュポットの速さは,
d(x − x1 )/dt になる.
CS を通して伝わる力は, F0 と R の並列を通して伝わる力に等しい.
dx dx1
x1
−
=R
+ F0
CS
dt
dt
17
(2.36)
図 2.6: 骨格筋のモデル
第二原理, つまり 2 つの並列からの力の和が F になることから,
F =
x1
x
+
CS
CP
(2.37)
式 (2.36) と式 (2.37) から x1 を消去する.
x1
CS
x1
dx1
dt
より, F , x, F0 の式をえる.
1
dF
+
F =
dt
RCS
= F−
x
CP
CS
x
CP
CS dx
dF
= CS
−
dt
CP dt
= CS F −
1
1
+
CS
CP
1
F0
dx
+
x+
dt
RCS CP
RCS
(2.38)
(2.39)
(2.40)
(2.41)
定常の等尺性収縮では (x = 0, dx/dt = 0, dF/dt = 0), F = F0 になる.
2.4
分布定数モデルと集中定数モデル
ここまで述べてきたモデルは, 集中定数モデル. 1 つの要素に集約されることを仮定している. 例
えば, 肺のモデルでは, 中央の気道の抵抗は, 単一の抵抗 RC になっている. 実際には, 異なる流体
の抵抗をもつ分岐がある. CL も単一の値である. 実際には, 肺の場所によって異なる. 空間的な分
布を, より現実的にする, 分布定数モデルがしばしば有用である. 空間 (長さあるいは体積) と時間
の偏微分方程式になる.
18
図 2.7: 無随神経のケーブル特性の集中定数系と分布定数系モデル
19
分布定数モデルは, たくさんの微小な集中定数のサブモデルのネットワークとみなすことができ
る. 無随神経線維の受動ケーブル特性のモデルを導く. 図 2.7 に示すように, 神経線維は, 直列につ
ながったサブユニットのネットワークとしてモデル化できる. サブユニットは, 回路要素 rx , rM ,
CM からなる. rx は, 断面積 1cm2 あたり, 長さ 1cm あたりの軸方法の抵抗 (単位は Ωcm) である.
rM と CM は, 神経膜表面の 1cm2 あたりの抵抗とキャパシタンスである. 神経線維が浸かっている
外液を電気的な接地と仮定する. このケーブルの任意の点の電位で流れる電流をどのように関連づ
けるのか.
分布定数モデルでは, 神経線維は, 半径 a の長い円筒状の形で考える. 微小な長さ δx を考える.
δx の間に, 細胞内電位が δV 上がって, 軸方向の電流が δI 増えたと仮定する (実際には, 電流が洩
れのある膜をとおして神経線維の外に流れるので, 神経線維に沿って電位は下がる). 長さ方向の V
の変化は負であることに注意して
Irx δx
(2.42)
πa2
式 (2.42) の右辺の負は, 電位の増加は電流の流れの方向と反対だから. 式 (2.42) の両辺を δx で割っ
δV = −
て, δx を 0 にする極限をとることによって,
Irx
∂V
=− 2
∂x
πa
(2.43)
膜電流 δI は, 細胞内電位 V に次のように関係する.
V
∂V
δI = −
2πaδx + cM 2πaδx
rM
∂t
(2.44)
式 (2.44) の両辺を δx で割って,
∂I
=−
∂x
式 (2.43) を x で偏微分し,
と式 (2.45) から,
V
∂V
2πa + cM 2πa
rM
∂t
∂I
∂ 2 V πa2
=− 2
∂x
∂x rx
2rx
∂2V
=
2
∂x
a
V
∂V
+
cM
∂t
rM
(2.45)
(2.46)
(2.47)
この式は, 一次元のケーブル方程式.
2.5
線形システムと重ねの理
線形システムは, 線形の常微分方程式あるいは偏微分方程式で特徴づけられる. 入力, 出力の変
数とその微分が 1 次であるとき, 微分方程式は線形である.
d2 y
dy
d2 x
dx
+
y
=
b
+ b0 x
+
a
+ b1
1
2
dt2
dt
dt2
dt
y とその微分, x とその微分に関する項のベキが全て 1 乗なので, 線形.
係数 a1 , a2 , b0 , b1 , b2 が定数であれば, 式 (2.48) は, 時不変 (time invaruant) である.
a2
(2.48)
一方, もしこれらの係数の 1 つ以上が, t の関数であったら, 例えば, a1 = t であったら, 式 (2.48)
はまだ線形ではあるが, 時変 (time varying) である.
しかし, これらの係数の 1 つ以上が, y あるいは x の関数であると, 式 (2.48) は, 非線形になる.
20
入力の経時変化を x = x1 (t) とするとき, 出力の経時変化が y = y1 (t) になったと考える.
a2
d2 y1
dy1
d2 x1
dx1
+
a
=
b
+ b1
+
y
+ b0 x1
1
1
2
dt2
dt
dt2
dt
(2.49)
入力の経時変化が別の x = x2 (t) とするとき, 出力の経時変化が y = y2 (t) になったと考える. 同
様に, 式 (2.50) を得る.
a2
d2 y2
dy2
d2 x2
dx2
+
y
+ b0 x2
+
a
=
b
+ b1
1
2
2
2
2
dt
dt
dt
dt
(2.50)
式 (2.49) と式 (2.50) を足す.
a2
d2 (y1 + y2 )
d(y1 + y2 )
d2 (x1 + x2 )
d(x1 + x2 )
+ (y1 + y2 ) = b2
+ b0 (x1 + x2 ) (2.51)
+ a1
+ b1
2
dt
dt
dt2
dt
2 つの異なる入力の和に対する線形システムの応答は, 個々の入力に対する応答の和である.
2 つ以上の入力の拡張して, 重ねの理として知られている. 線形システムの定義であり, システム
が線形であるか否かのテストにしばしば用いられる.
重ねの理は, 式 (2.48) の完全な解は, 2 つの部分にわけられることを示す.
y(t) = yc (t) + yp (t)
(2.52)
yc (t) は余関数 (同次方程式の解) で, yp (t) は特解と呼ばれる.
yc (t) は, 式 (2.48) の入力が 0 のときの線形システムの応答である.
a2
d2 yc
dyc
+ a1
+ yc = 0
dt2
dt
(2.53)
yc (t) は, y(t) の内, x(t) に独立な成分を表す. しかし, yc (t) は, y およびそれらの微分で与えられる
初期値に依存する.
yp (t) は, 入力 x(t) を, 時間の特別な関数にしたときの同じ線形システムの応答である.
a2
d2 yp
dyp
d2 x
dx
+
y
+ b0 x
+
a
=
b
+ b1
1
p
2
dt2
dt
dt2
dt
(2.54)
yp (t) は, 式 (2.54) の係数と入力 x(t) の経時変化にのみ依存する.
もし, 線形システムが安定であれば, yc (t) で特徴付けられる y(t) の成分は, 減衰して零になり, 全
体の応答は yp (t) が支配的になる.
x(t) がずっと持続しているときに, yc (t) はしばしば過渡応答と呼ばれ, yp (t) は定常状態の応答
として知られている.
2.6
ラプラス変換と伝達関数
線形システムの入出力関係は, 微分方程式で表すことができる. 数学的記述とブロック線図を結
び付ける.
ラプラス変換は以下のとおり.
L[y(t)]
≡
Y (s) ≡
L−1 [Y (s)]
≡
1
2πj
21
∞
0
σ+j∞
σ−j∞
y(t)e−st dt
(2.55)
Y (s)est ds
(2.56)
ラプラスの変数 s は複素数である (s = σ + jω).
式 (2.55) を y(t) に適用すると, 複素数の s 空間の関数 Y (s) を得る.
y(t) の時間微分, dy/dt に適用すると, 次のようになる.
∞
∞
dy(t)
dy(t) −st
L
e dt = s
y(t)e−st dt + [y(t)e−st ]∞
≡
0
dt
dt
0
0
よって,
(2.57)
dy
L
= sY (s) − y(0)
dt
同じ変換を d2 y/dt2 に適用する.
2 dy
d y
2
L
= s Y (s) − sy(0) −
2
dt
dt t=0
y(t) の積分に適用する. ただし, t < 0 のとき y(t) = 0 とする.
t
Y (s)
L
y(t)dt =
s
0
(2.58)
(2.59)
(2.60)
式 (2.48) にラプラス変換を適用する. ただし, t = 0 のとき, x, y およびそれらの微分は 0 と仮
定する.
a2 s2 Y (s) + a1 sY (s) + Y (s) = b2 s2 X(s) + b1 sX(s) + b0 X(s)
(2.61)
初期値が 0 でなくても, 線形システムの動的な関数の性質には影響しない.
式 (2.61) を次のように変形する.
b 2 s2 + b 1 s + b 0
Y (s)
=
X(s)
a2 s2 + a 1 s + 1
(2.62)
どんな入力 x(t) についても, そのラプラス変換 X(s) を計算することができて, それを式 (2.62)
の右辺に乗ずることによって, システムの応答のラプラス変換 Y (s) を得ることができる. Y (s) を
逆変換することで, 時間領域の応答 y(t) を復元することができる.
ラプラス変換は, 微分方程式を代数方程式に変換する. 代数方程式は解きやすい.
Y (s) と X(s) の比は, 伝達関数と呼ばれる. ブロック線図でかける.
逆ラプラス変換する前に, 伝達関数を部分分数に展開する. 例えば, 式 (2.62) で, b2 = 0, b1 = 0.5,
b0 = 1.25, a2 = 0.25, a1 = 1.25 とする. 入力 x(t) を単位ステップ関数とすると, X(s) = 1/s.
1
0.5s + 1.25
1
1
1
Y (s) =
=
−
−
+ 0.25
(2.63)
s(s + 1)(0.25s + 1)
s s+1
s s+4
時間領域の出力 y(t) は,
y(t) = (1 − e−t ) + 0.25(1 − e−4t )
2.7
(2.64)
インパルス応答と畳み込み積分
未知の伝達関数 H(s) の線形システムを考える. Y (s) = H(s)X(s) なので, ラプラス変換が 1 に
なる入力 x(t) をみつけることができれば, Y (s) = H(s) なので, 出力のラプラス変換を未知のシス
テムの伝達関数が明らかになる.
22
表 2.1: ラプラス変換表
x(t)
X(s)
δ(t)
u(t)
1
e−at u(t)
t −at
u(t)
k! e
k
tu(t)
tk
k! u(t)
sin(ωt)u(t)
cos(ωt)u(t)
e−at sin(ωt)u(t)
e−at cos(ωt)u(t)
e
−at
cos(ωt + θ)u(t)
t sin(ωt)u(t)
t cos(ωt)u(t)
1
s
s
s+a
1
(s+a)k+1
1
s2
1
sk+1
ω
s2 +ω 2
s
s2 +ω 2
ω
(s+a)2 +ω 2
s+a
(s+a)2 +ω 2
(s+a) cos(θ)−ω sin(θ)
(s+a)2 +ω 2
2ωs
(s2 +ω 2 )2
2
s −ω 2
(s2 +ω 2 )2
X(s) = 1 になる時間領域の関数は, Dirac の δ(t) である. 無限大の振幅と無限小のパルス幅をも
つ矩形パルス.
X(s) = 1 のとき, Y (s) = H(s) になり, 逆ラプラス変換して, y(t) = h(t). インパルスが入力さ
れたときの応答, つまりインパルス応答 h(t) は, 伝達関数 H(s) の逆ラプラス変換である.
H(s) が既知の伝達関数で, X(s) が任意の入力 x(t) のラプラス変換である場合を考える.
y(t) =
=
=
=
L−1 [H(s)X(s)]
σj +∞
1
H(s)X(s)est ds
2πj σ−j∞
σ+j∞ ∞
1
h(τ )e−sτ dτ · X(s)est ds
2πj σ−j∞ 0
∞
σ+j∞
1
X(s)es(t−τ ) ds · h(τ )dτ
2πj
0
σ−j∞
(2.65)
(2.66)
(2.67)
(2.68)
内側の積分について考える.
逆ラプラス変換の定義から,
1
2πj
σ+j∞
σ−j∞
X(s)es(t−τ ) ds = x(t − τ )
つまり, 式 (2.65) は, 次のようになる.
y(t) = L−1 [H(s)X(s)] =
0
∞
h(τ )x(t − τ )dt
(2.69)
(2.70)
式 (2.70) は, 畳み込み積分. 線形システムにおいて, もっとも重要でもっとも基本的な演算.
• 時間領域の乗算は, s 領域の畳み込み積分. 時間領域の畳み込み積分は, s 領域の乗算.
• インパルス応答 h(t) は, システムの動的な特性を与える. h(t) が既知であれば, 任意の入力
x(t) に対する応答 y(t) を, x(t) と h(t) の畳み込み積分によって導くことができる.
23
• 線形システムのダイナミクスは, 伝達関数 H(s) で完全に表される. 伝達関数 H(s) が分かれ
ば, 任意の入力に対するシステムの応答のラプラス変換 Y (s) を, 入力のラプラス変換 X(s)
と伝達関数 H(s) の積で求めることができる.
x(t) = u(t) のとき, 式 (2.70) は, 次のようになる.
t
y(t) = g(t) =
h(τ )dτ
(2.71)
0
ここで, g(t) はステップ応答.
2.8
状態空間による解析
生体は複雑なので, 高次の微分方程式が必要になる. 複雑な微分方程式の解析解は, 一般的には
得られないし, 数値解は不安定さの問題に悩まされる. ラプラス変換は便利な方法だが, 初期値が
零ではない場合には, 複雑になることがある. このような条件では, 状態空間でのモデリングが有
効である. 状態空間モデルは, 時変システムや非線形システムへの拡張が容易である.
N 次の微分方程式で表せる線形システムを考える.
dN −1
dy
dN y
+ a0 y = b0 x(t)
+
a
+ · · · + a1
N
−1
dtN
dtN −1
dt
(2.72)
N 個の状態変数 z1 (t), z2 (t), · · · , zN (t) を定義する.
z1 (t)
=
z2 (t)
=
y(t)
dz1 (t)
dy(t)
=
dt
dt
(2.73)
(2.74)
..
.
(2.75)
N −1
zN (t)
d
=
y(t)
dzN −1 (t)
=
dtN −1
dt
(2.76)
式 (2.72) は, 次のように書ける.
dzn
= −a0 z1 − a1 z2 − · · · − aN −1 zN + b0 x
dt
(2.77)
行列の微分方程式でかける.
dz(t)
= F z(t) + Gx(t)
dt
(2.78)
ここで
⎡
⎢
⎢
⎢
⎢
z(t) = ⎢
⎢
⎢
⎣
z1 (t)
z2 (t)
..
.
zN −1 (t)
zN (t)
⎤
⎡
⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎥
⎥,F = ⎢
⎢
⎥
⎢
⎥
⎣
⎦
0
0
..
.
0
1
0
..
.
0
−a0
−a1
0 ·
1 0
.. ..
. .
0 ·
·
·
·
·
..
.
0
0
0
..
.
1
−aN −2
−aN −1
⎤
⎡
⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎥
⎥,G = ⎢
⎢
⎥
⎢
⎥
⎣
⎦
0
0
..
.
0
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(2.79)
b0
高次の微分方程式を, 1 次の微分方程式 N 個に変換する.
状態ベクトル z(t) と出力 y(t) は, 次の式で関係付けられる.
y(t) = Cz(t) + Dx(t)
24
(2.80)
ここの例では,
C =[ 1 0
· · · 0 0 ], D = 0
(2.81)
式 (2.78) を状態方程式, 式 (2.80) を観測方程式と呼ぶ.
線形システムでは, 伝達関数から状態方程式, 状態方程式から伝達関数への変換は容易である.
式 (2.62) を例に計算する. 中間の変数 U (s) を定義する.
H(s) =
1
b 2 s2 + b 1 s + b 0
U (s) Y (s)
=
2
X(s) U (s)
a2 s + a1 s + 1
1
(2.82)
a2 s2 U (s) + a1 sU (s) + U (s) = X(s)
(2.83)
積の第 1 項から
逆ラプラス変換して,
d2 u(t)
du(t)
+ u(t) = x(t)
+ a1
2
dt
dt
状態変数 z1 (t) および z2 (t) を次のように定義する.
a2
(2.84)
du(t)
dz1 (t)
=
dt
dt
(2.85)
dz2 (t)
a1
1
x(t)
= − z2 (t) − z1 (t) +
dt
a2
a2
a2
(2.86)
z1 (t) = u(t), z2 (t) =
式 (2.84) は, 次のように書ける.
式 (2.85) と式 (2.86) を合わせて, 次の行列の状態方程式を得る.
⎡
⎤ ⎡
⎤
⎤
⎡
dz1 (t)
0
1
0
(t)
z
1
⎢ dt ⎥ ⎣
+ ⎣ 1 ⎦ x(t)
1
a1 ⎦
⎣ dz (t) ⎦ =
2
−
−
z2 (t)
a2
a2
a2
dt
(2.87)
H(s) の残りの項は次のようになる.
Y (s) = b2 s2 U (s) + b1 sU (s) + b0 U (s)
(2.88)
逆ラプラス変換して, 次の式を得る.
y(t) = b2
d2 u(t)
du(t)
+ b0 u(t)
+ b1
dt2
dt
(2.89)
式 (2.85) および式 (2.86) を用いて,
z1
a1
b2
b2
+ x(t)
y(t) =
b0 −
b1 − b2
a2
a2
a
z2
2
2.9
(2.90)
計算機による解析とシミュレーション —MATLAB と SIMLINK—
ラプラス変換や状態空間によるモデル化は, 線形システムの数学的特性を簡単化するが, 現実の
動的な振舞を与えるモデルは, 複雑過ぎて解析的には扱えない. このように複雑な場合には, ブロッ
ク線図で表されたシステムを, 計算機モデルに変換し, 数値的に解くことである. 従来, 微分方程式
25
は, 基本的なプログラミング言語である C や Fortran で開発されたプログラムで解かれた. 現在,
沢山尾のソフトウエアを利用可能で, もっと簡単にモデルシミュレーションを行うことができるも
のがある. それらの 1 つに SIMLINK がある. SIMLINK は, グラフィカルな環境を提供し, ユー
ザは, ブロック線図を数学関数のブロックのネットワークに容易に変換することができる. コマン
ドベースの会話的な環境である MATLAB とともに実行する. 少なくとも, 学生版の MATLAB と
SIMLINK を利用できること, MATLAB の基本的な機能を知っていることを前提にする.
人工呼吸器のピーク圧が定められたときに, 患者にどれだけの体積を与えられるか調べたいと仮
定する. この問題の解は, 患者の肺の機構に関する知識を必要とする. この患者は, 比較的正常な
機能をもっており, 肺のパラメータを以下のとおりに仮定する. Rc = 1 cm H2 OsL−1 , Rp = 0.5 cm
H2 OsL−1 , CL = 0.2 L H2 O−1 , Cw = 0.2 L H2 O−1 , Cs = 0.005 L H2 O−1 . 2 つのアプローチを
とる.
1 つは, 直接的な方法で, 全体の伝達関数を導くことであり, SIMLINK の 1 つのブロックとして
用いる. 全体の空気の流れ Q と気道の入口の圧力 Pao に関する全体の微分方程式は, キルヒホフの
法則で導かれ, 式 (2.34) が得られる. 上記のパラメータを代入して, ラプラス変換すると, 次の式が
得られる.
s2 + 420s
Q(s)
= 2
Pao (s)
s + 620s + 4000
(2.91)
1. 上記のモデルを実装するために, MATLAB コマンドウインドウで SIMULINK を実行する
(MATLAB プロンプトで simulink と入力する). SIMULINK のライブラリウインドウが表示
される.
2. File メニューから New の model を選ぶ. untitled の別のウインドウが表示される. これが
我々のモデルを作成する作業ウインドウである.
3. Continuous を開いて, Transfer fnc を選び, 作業用ウインドウへドラッグする.
4. このブロックをダブルクリックして, 伝達関数のパラメータを入力する. 分子は [1 420 0], 分
母は [1 620 4000].
5. Pao を作成する. Source を開いて, sine wave を選び, 作業用ウインドウへドラッグする.
6. ダブルクリックしてパラメータを入力する. 振幅 2.5 cm H2 O, 各週は数 1.57 rad s−1 . これ
は, 0.25 Hz(一分間に 15 回) に相当する.
7. 出力を伝達関数の入力につなぐ.
8. 出力 Q を観測するために, Sinks から Scope を選び, 作業用ウインドウにドラッグする.
9. Scope を Pao に接続する.
10. Scope をダブルクリックして, 横軸・縦軸の範囲などを設定する.
11. 別の Scope を作業用ウインドウにドラッグして, Q と接続する. 軸の範囲などの設定を同様
に行う.
12. 患者の肺に入る体積を求めるために, Conitinuous から integrator を選び, 作業用ウインドウ
にドラッグする.
13. Integrator に Q を接続し, 積分した結果を別の Scope に接続する.
26
14. シミュレーション結果をファイルに出力するために, Sinks から To File を選び, 作業用ウイ
ンドウにドラッグする. ダブルクリックして, ファイル名を respm1.mat にする.
15. Pao , Q, 体積をファイルに出力するので, 3 つをまとめる Mux を作業用ウインドウにドラッ
グする. Mux をダブルクリックして, 入力数を 3 に変更する.
16. Pao , Q, 体積を Mux の入力に接続し, Mux の出力を To File に接続する.
次に, シミュレーションを実行する.
1. Simulation メニューから Parameters を選ぶ. Solver で, Fixed-step, 0.001 に設定する／
2. Simulation → Start を実行する.
2 つ目は, 伝達関数の形に変換して SIMULINK に実装するのではなく, 内部の状態も分かるよう
に実装する.
ブロック線図の要素が, 式 (2.25) と式 (2.26) に直接的に対応する.
Paw = RP QA +
1
1
+
CL
CW
QA dt
PAO
1
Cs
(Q − QA )dt
1
= RC Q +
(Q − QA )dt
CS
=
1/CL , 1/Cw , Rp , Cs を含む 2 重ループは, 式 (2.25) を表す. この実装では, 個々のパラメータが
変化するときに, 何が起こるかを調べることが容易である. しかし, 代数ループができるので, 特別
な対策が必要になることがある.
27