復習 ) 時系列のモデリング ~a. 離散時間モデル ~ y k + a 1 z 1 y k + + a na z n ay k = b 0 u k + b 1 z 1 u k + + b nb z n bu k y k = G z 1 u k = B(z 1 ) A(z 1 u k ) ARMA モデル A z 1 B z 1 = 1 + a 1 z 1 + + a na z n a = b 0 + b 1 z 1 + + b nb z n b 連続時間モデルと離散時間モデルは z 変換によって関連付けられる. 計算機へのプログラミングや実装では, 未来のデータを利用することはできない. 過去のデータ系列を明記するために z 1 演算子を利用する. 伝達関数の形, 誤差の伝達関数の形によって, 多様なモデルがある. e.g. AR, MA, ARIMA, CARMA, etc. 1
復習 ) 時系列のモデリング ~a. 離散時間モデル ~ ARMA モデルの様々な表現形式 y k = B(z 1 ) A(z 1 u k ) スカラ関数 = b 0 + b 1 z 1 + + b nb z n b 1 + a 1 z 1 + + a na z n a = β 0 + β 1z 1 + + β nb z nb 1 + a 1 z 1 + + a na z n a n y k = a i=1 n a i y k i + b j=0 b j u k j 1 可観測正準形 x 1 k + 1 x k + 1 x n k + 1 = 0 0 a na 1 0 0 a na 1 0 1 0 a na 0 0 1 a 1 y k = 0 0 1 x 1 k x k x n k x 1 k x k x n k x(k + 1) A x(k) b c T x(k) + β 0 u(k) 直達項 + β 1 β β n u(k)
復習 ) 時系列のモデリング ~a. 離散時間モデル ~ 実用に際しては, ホールドの伝達関数を考慮する必要がある. ZOH(zero-order hold) の場合 y k = H z 1 G z 1 u k = z 1 B(z 1 ) A(z 1 u k ) ZOH があれば直達項は存在しない y(k) = b 0z 1 + b 1 z + + b nb z n b 1 1 + a 1 z 1 + + a na z n a = β 0z 1 + + β nb z n b 1 1 + a 1 z 1 + + a na z n a 3
3. システムのモデリング 4
3.1 信号とシステム モデリング (modeling) 入力信号 u と出力信号 y が与えられたときシステムの特性 G を求める. 解析 (analysis) 入力信号 u とシステムの特性 G が与えられたとき出力信号 y の性質を求める. 設計 (design) システムの特性 G と出力信号 y の目標値 r が与えられたとき入力信号 u を求める. これれらの語句の正確な違いを理解するのは重要 5
3. 制御のためのモデリング モデルベースト制御 (MBC: model-based control) 現代制御 ロバスト制御 モデル予測制御 etc. モデルフリー制御 (MFC: model-free control) PID 制御 ファジィ制御 ニューロ制御 etc. 近年, 産業界では MBC の重要性が認識され始めている. 計算機の発展と低価格化によって,MBC の実装が容易になっている. 6
3. 制御のためのモデリング ホワイトボックスモデリング (white-box modeling) 第一原理モデリング (first principle modeling) や物理モデリングとも呼ばれる. 対象を支配する第一原理 ( 運動方程式, 回路方程式, 電磁界方程式, 保存則, 化学反応式などの科学原理 ) に基づいてモデリングを行う方法. 利点 対象の挙動が忠実に再現できる. 計算機上での解析や予測に利用するシミュレータの構築に役立つ. 問題点 一般に, 非線形 偏微分方程式で記述され複雑である. 制御系設計にそのまま利用することは難しい. 実験を行わなければ正確な値が分からない ( 実験が困難な場合も多い ) パラメータが存在する. 7
3. 制御のためのモデリング ブラックボックスモデリング (black-box modeling) システム同定 (system identification) とも呼ばれる. 対象をブラックボックスと見なして, その入出力データから統計的な手法でモデリングを行う手法. 線形システム同定は論理体系が完備されている. ビッグデータからモデル化に有用な情報をいかに抽出するかは, 重要な課題の一つとして注目されている. 利点 複雑なシステムに対しても, 実験データから比較的簡潔なモデルを得ることができる. 様々なシステム同定法が提案されており, それを実行するソフトウェアも MATLAB に用意されている. 問題点実験的モデリング手法であるので, モノが無いとモデリングできない. 8
3. 制御のためのモデリンググレーボックスモデリング (grey-box modeling) 部分的な物理情報が利用できる場合のモデリング法. 利点 最も現実的なモデリング法であり, 実際の制御の現場で用いられるモデリング法のほとんどは, この範疇に含まれる. 問題点 対象や実験環境に大きく依存するため, グレーボックスモデリングの一般的な方法論を構築するのは難しい. 9
3.5 制御のためのモデリングのポイント 実システム 詳細モデル ノミナルモデル 実システムをできるだけ忠実に再現する詳細モデルを構築しようとする立場 ( 第一原理モデリング ) ノミナルモデルと実システムの差が モデルの不確かさ 制御用ノミナルモデルはできるだけ簡単なものが望ましい 1. どれだけ第一原理モデルを実際の制御対象に近づけられるか?. 制御用ノミナルモデルが第一原理モデルの重要な部分をよく近似しているか? 3. 近似しきれなかった部分を定量的に評価できるか? 不確かさのモデリング と呼ばれる研究分野 10
4. スカラ変数とベクトル変数に対する最小二乗法 11
4.1 最小二乗推定法 ( スカラの場合 ) 最小二乗推定値の導出 確率変数 : x(k) 観測雑音 : w(k) 観測 : y k = cx k + w(k) 既知確率変数の平均値 : E x(k) = x 確率変数の分散 : E x k x = σ x 観測雑音の平均値 : E w(k) = w 観測雑音の分散 : E w k w = σ w 推定問題 : x k = αy k + β (α, β) を求める. 条件 : e k = x k x k の平均二乗誤差の最小化 1
4.1 最小二乗推定法 ( スカラの場合 ) 1 1 次モーメントに関する条件 推定誤差の平均値を最小化する (α, β) の関係性は? どのようなものだろう? E e(k) = E x(k) x = E x k αy(k) β = E x k α cx k + w k β = 1 αc = 0 x α w β β = 1 αc x α w = x α c x + w 13
4.1 最小二乗推定法 ( スカラの場合 ) 次モーメントに関する条件 推定誤差の分散を最小化する (α, β) の関係性は? どのようなものだろう? E e k E[e(k)] = E x(k) αy(k) β 1 αc x α w β = E x(k) α cx(k) + w(k) 1 αc x + α w = 1 αc E x(k) x 1 αc αe x k x w k w +α E w(k) w 0 = 1 αc σ x + α σ w さらに α についての平方完成により = c σ x + σ w α cσ x α + σ x cσ x = c σ x + σ w α c σ x + σ w cσ w = c σ x + σ w α σ x + c σ w cσ x 4 + σ x c σ x + σ w 1 + σ x + c σ w 14
4.1 最小二乗推定法 ( スカラの場合 ) cσ w E e k E[e(k)] = c σ x + σ w α σ x + c σ w 1 + σ x + c σ w = c σ x + σ w α cσ w σ + σ σ これより,α = cσ w σ のとき, 推定誤差分散は最小値 E e k E[e(k)] = σ をとる. また, x k の推定値は, x k = αy k + β = cσ w σ y k + x cσ w σ c x + w = x + cσ w σ y k c x + w で与えられる. 15
4.1 最小二乗推定法 ( スカラの場合 ) 直交性の原理 最小二乗推定値 x(k) と推定誤差 e(k) の相関関数を計算してみる E x k e k = E x + cσ σ w y k c x + w e k = x cσ σ w c x + w E e(k) + cσ σ w E[y(k)e(k)] = cσ σ w E[y(k) x(k) x(k) ] ここで E e(k) = 0 を利用した 次に式の中の E[y(k) x(k) x(k) ] を計算する 16
4.1 最小二乗推定法 ( スカラの場合 ) E y(k) x(k) x(k) = E cx(k) + w(k) x(k) x + cσ σ w y(k) c x + w = cx k + w k 1 cσ σ w x k x cσ σ w w k w = c 1 cσ E x k x k x σ w cσ σ w E w k w k w となる. ここで, E x(k) x(k) x = E x(k) x x(k) x + x x(k) x = σ x E w(k) w(k) w = E w(k) w w(k) w + w w(k) w = σ w を代入すると, 17
4.1 最小二乗推定法 ( スカラの場合 ) E y k x k x k = c 1 cσ σ w σ x cσ σ w σ w = c 1 σ c σ w + 1 σ x + σ σ x σ x cσ となり, が成り立つことが導かれた = c 1 σ 1 σ + σ σ x σ x cσ = 0 E e(k)y(k) = 0 この式は, 最小二乗推定誤差と観測は無相関である ことを示す. これは直交性の原理とよばれる 18
4. 最小二乗推定法 ( 多変数の場合 ) 対象 y k = Cx k + w(k) 確率変数 : x(k) R n 観測雑音 : w(k) R n 観測ベクトル : y(k) R n 観測行列 : C R n n 既知 x(k) とw(k) は無相関確率変数の平均値 : E x(k) = x 確率変数の分散 : E (x k x) x k x T = Σ x 観測雑音の平均値 : E w(k) = w 観測雑音の分散 : E (w k w) w k w T = Σ w (4.39) (4.40) 19
4. 最小二乗推定法 ( 多変数の場合 ) 出力の平均ベクトル E y = E Cx + w = C x + w 共分散行列 E (y E y )(y E y ) T = E (Cx + w C x w)(cx + w C x w) T = E C x x + (w w) {C x x + (w w)} T = CΣ x C T + Σ w 推定誤差ベクトル e k = x k x k 推定問題 x k = Fy k + d (F R n n, d R n ) を求める. (4.4) 0
4. 最小二乗推定法 ( 多変数の場合 ) 1 1 次モーメントに関する条件? 推定誤差の平均値を最小化する d はどのようなものだろう? E e = E x Fy d = x F C x + w d = 0 d = x F C x + w = I FC x F w (4.45) 次モーメントに関する条件? 推定誤差の分散を最小化する F はどのようなものだろう? P E ee T 誤差共分散行列 P = E[(x Fy d)(x Fy d) T ] = E[{x F Cx + w d}{x F Cx + w d} T ] (4.47) 1
4. 最小二乗推定法 ( 多変数の場合 ) x F Cx + w d = I FC x Fw d = I FC x Fw I FC x + F w = I FC (x x) F(w w) Eq. (4.45) を利用 (4.48) Eq. (4.48) を Eq. (4.47) に代入 P = E[{ I FC (x x) F(w w)}{ I FC x x F w w } T ] = I FC E x x x x T I FC T + FE[ w w w w T ]F T = I FC Σ x I FC T + FΣ w F T Eq. 4.39 (4.40) を利用 = F CΣ x C T + Σ w F T FCΣ x Σ x C T F T + Σ x A B (4.50) P = FAF T FB B T F T + Σ x (4.51)
4. 最小二乗推定法 ( 多変数の場合 ) P = F B T A 1 A(F B T A 1 ) T +Σ x B T A 1 B F = B T A 1 のとき P は最小 F = Σ x C T (CΣ x C T + Σ w ) 1 (4.55) P = Σ x B T A 1 B (4.54) Eq. 4.45,(4.55) をEq. 4.4 に代入 x k = Fy + I FC x F w = x + F{y (C x + w)} = x + Σ x C T (CΣ x C T + Σ w ) 1 {y (C x + w)} 最小二乗推定ベクトル 3
4. 最小二乗推定法 ( 多変数の場合 ) P = Σ x Σ x C T (Σ w + CΣ x C T ) 1 CΣ x Eq. 4.54 (4.50) を利用 = (Σ x 1 + C T Σ w 1 C) 1 (4.60) Eq. 4.55 を変形 F = PC T Σ w 1 (4.6) 最小二乗推定 x k = x + F{y(k) (C x + w)} F = Σ x C T (CΣ x C T + Σ w ) 1 P = Σ x Σ x C T (Σ w + CΣ x C T ) 1 CΣ x 4