はじめに
状態方程式によるモデル化でシステムの振る舞いを数学的に表現した状態方程式にまとめました。
\[\frac{dx(t)}{dt}=Ax(t)+Bu(t)\]この状態方程式をデジタル制御システムで扱えるように離散化した逐次式
\[x(t+\Delta t)=P x(t)+Q u(t) \]に変換します。
システムを離散化した逐次式で表現する理由
システムの制御ロジックは一般的にデジタル制御システム（マイコンを搭載したコントロールユニット）に実装するため、下図のような離散時間制御になります。また、離散化した逐次式にすると、Excelでも時間応答を簡単に計算でき、対象となるシステムの振る舞いを確認（シミュレーション）できるようになります。
組み込み制御のイメージ
状態方程式を離散化した逐次式に変換
状態方程式を離散化した逐次式に変換する手法を説明しますが、途中に出てくる微分方程式の一般解やPade近似などは参考知識なので、具体的な手法は最後の PとQの計算式 になります。
状態方程式の解を求める
微分式のままではシステムの振る舞いを具体的に計算することが難しいので、状態方程式の解からシステムの時間応答を計算できるようにします。
状態方程式は連続時間の微分方程式
\[\frac{d x(t)}{d t}=A x(t)+B u(t) \\\]微分方程式の一般解は
\[\boldsymbol{x}(t)=e^{A t} x_{0}+\int_{0}^{t} e^{A(t-\tau)} B u(\tau) d \tau \\\]
離散式の導出
この一般解を使って離散式を導き出します。
時刻tにおける値を𝒙(𝑘) としたとき、Δt 後の t+Δt の状態 𝒙(𝑘+1) は
\[\boldsymbol{x}(k+1)=e^{A \Delta t} \boldsymbol{x}(k)+A^{-1}\left(e^{A \Delta t}-I\right) B \boldsymbol{u}(k) \]ここで \(P=e^{A \Delta t} \\\) \(Q=A^{-1}\left(e^{A \Delta t}-I\right) B \) とすると
\[x(t+\Delta t)=P x(t)+Q u(t) \]
になります。
線形システムで、状態方程式のA,Bで表すシステムのパラメータが一定なら、離散時間 Δt を定めれば、PとQを計算することで離散式に変換できます。
\(e^{A \Delta t} \)の計算方法
一般にAは行列なので\(e^{A \Delta t} \)の直接計算は難しいため、このブログでは Pade近似 を使います。
Pade近似式でPとQを計算すれば状態方程式は逐次式に変換することができます。
Pade近似 \[e^{x}=\frac{1+x / 2}{1-x / 2}\]を適用すると、P Q は以下の式になります。
\[P=e^{A \Delta t}=\left(I+\frac{\Delta t}{2} A\right) \quad\left(I-\frac{\Delta t}{2} A\right)^{-1}=I+A \Delta t\left(I-\frac{\Delta t}{2} A\right)^{-1} \]\[Q=A^{-1}\left(e^{A \Delta t}-I\right) B=\Delta t\left(I-\frac{\Delta t}{2} A\right)^{-1}B\]
Excelには行列演算の関数 MINVERSE や MMULT があるので、PとQはExcelで容易に計算できます。計算方法は次のCR回路モデルの時間応答の確認で具体的なモデルを例に説明していきます。
最後までお読みいただき、ありがとうございました。