微分方程式

微分方程式の基礎

微分方程式はある関数とそれを微分した導関数の関係式であり，関数の特定の変数に対する変化を記述することができる．まず，1階線形微分方程式を例として見てみよう．

\[ \frac{dx(t)}{dt}=a_c x(t)+b_c u(t) \]

状態変数 \(x(t)\)は，時間\(t\)に対する関数である．

添え字の\(c\)は連続 (continuous) を意味するが，これは後で離散化する際に区別するためである．この方程式においては\(b_c=0\)の場合を同次方程式, \(b_c\neq 0\)の場合を非同次方程式という．

微分方程式の解

微分方程式を解くとは\(x(t)\)のような関数の具体的な式を求めることである．上式の解は

\[ x(t)=e^{a_c t}x(0)+\int_0^t e^{a_c (t-\tau)}b_c u(\tau) d\tau \]

として与えられる．微分方程式を解く手法は様々で，それぞれの方程式について適切な手法を選択する．本書ではLaplace変換を多用するが，細かい説明は付録にて行う．

連立線形微分方程式

\(n\)個の微分方程式

連立線形微分方程式という．これをベクトル，行列を用いて

時不変 (time-invariant) の定数行列を\(\mathbf{A}_c \in \mathbb{R}^{n\times n}, \mathbf{B}_c \in \mathbb{R}^{n\times m}\), 状態ベクトルを\(\mathbf{x}(t)\in\mathbb{R}^n\), 入力ベクトルを\(\mathbf{u}(t)\in\mathbb{R}^m\)とする．

\[ \frac{d\mathbf{x}(t)}{dt} = \mathbf{A}_c\mathbf{x}(t) + \mathbf{B}_c\mathbf{u}(t) \]

解は

\[ \mathbf{x}(t)=e^{t\mathbf{A}_c}\mathbf{x}(0)+\int_0^t e^{(t-\tau)\mathbf{A}_c}\mathbf{B}_c\mathbf{u}(\tau) d\tau \]

連続時間モデルから離散時間モデルへの変換

離散化方法1: 解析解を用いた方法

1次元の場合

区間\([t, t+\Delta t]\)において入力\(u(t)\)が一定であると仮定すると，

\[\begin{split} \begin{aligned} x(t+\Delta t)&= \underbrace{e^{a_c \Delta t}}_{=: a_d}\mathbf{x}(t)+\underbrace{\left[\int_t^{t+\Delta t} e^{a_c(t+\Delta t-\tau)} d\tau\right] b_c}_{=: b_d}u(t)\\ &=a_d x(t)+b_d u(t)\\ \end{aligned} \end{split}\]

n次元の場合

区間\([t, t+\Delta t]\)において入力\(\mathbf{u}(t)\)が一定であると仮定すると，

\[\begin{split} \begin{aligned} \mathbf{x}(t+\Delta t)&=e^{(t+\Delta t)\mathbf{A}_c}\mathbf{x}(0)+\int_0^{t+\Delta t} e^{(t+\Delta t-\tau)\mathbf{A}_c}\mathbf{B}_c\mathbf{u}(\tau) d\tau\\ &=e^{\Delta t\mathbf{A}_c}e^{t\mathbf{A}_c}\mathbf{x}(0)+e^{\Delta t\mathbf{A}_c}\int_0^{t} e^{(t-\tau)\mathbf{A}_c}\mathbf{B}_c\mathbf{u}(\tau) d\tau + \int_t^{t+\Delta t} e^{(t+\Delta t-\tau)\mathbf{A}_c}\mathbf{B}_c\mathbf{u}(\tau) d\tau\\ &\approx \underbrace{e^{\Delta t\mathbf{A}_c}}_{=: \mathbf{A}_d}\mathbf{x}(t)+\underbrace{\left[\int_t^{t+\Delta t} e^{(t+\Delta t-\tau)\mathbf{A}_c} d\tau\right] \mathbf{B}_c}_{=: \mathbf{B}_d}\mathbf{u}(t)\\ &=\mathbf{A}_d\mathbf{x}(t)+\mathbf{B}_d\mathbf{u}(t)\\ \end{aligned} \end{split}\]

離散化した場合の

\[ \mathbf{x}_{t+1} = \mathbf{A}_c\mathbf{x}_t + \mathbf{B}_c\mathbf{u}_t \]

状態遷移方程式 (dynamics equations) とも呼ぶ．

離散化方法2: 微分方程式の数値解法

Euler法

Euler法は\(\dfrac{dx}{dt}=f(x, t)\)において，\(x_{n+1}=x_t+\Delta t f(x_n, t_n)\)とする手法である．

\[\begin{split} \begin{aligned} x(t+\Delta t)&=x(t) + \left[a_c x(t)+b_c u(t) \right]\Delta t\\ &=(1+a_c \Delta t)x(t) + b_c\Delta t u(t ) \end{aligned} \end{split}\]

ここで，解析解を用いる方法とEuler法の離散化係数の比較をしよう．

\(\Delta t=0\)でTaylor展開により1次近似すると

\[ e^{a \Delta t} \approx 1 + a\Delta t \]

\(a_c\neq 0\)の場合，

\[\begin{split} \begin{aligned} \int_t^{t+\Delta t} e^{a_c(t+\Delta t-\tau)} d\tau&=\frac{1}{a_c}(e^{a_c \Delta t}-1)\\ &\approx \frac{1}{a_c}\cdot a_c \Delta t=\Delta t \end{aligned} \end{split}\]

Runge-Kutta法

その他のsolver

adaptiveな方法など．JuliaであればDifferentialEquations.jlなどで実装されているsolverを用いる方が効率的である．

本書では主にEuler法を用いて実装を行う．Euler法は精度が低い手法であるという欠点があるものの，実装が簡便で可読性が高いことや，本書で扱うモデルに関してはEuler法でも定性的に同様の結果が再現できることなどが採用する理由である．