制御系CAD

制御を意識したモノつくりを目指して

１次系のLQ制御

投稿日時: 2021年10月8日投稿者: cacsd

１次系のLQ制御…Homework

[1] ある平衡状態周りの挙動が状態方程式

$\displaystyle{(1)\quad \dot{x}(t)=ax(t)+bu(t) }$

で表される制御対象に対して、状態フィードバック（state feedback）

$\displaystyle{(2)\quad \boxed{u(t)=-fx(t)} }$

を考えます。ここで、定数 $f$ はゲインと呼ばれ、 $a-bf<0$ を満たすものを選びます。(2)を(1)に代入して、次式を得ます。

$\displaystyle{(3)\quad \boxed{\dot{x}(t)=(a-bf)x(t)\quad (a-bf<0)} }$

これは制御対象とコントローラ（状態フィードバック）が閉路を構成することから（図1参照）、閉ループ系（closed-loop system）と呼ばれます。ここで、物理的に興味深いのは、制御対象が漸近安定でなくても（ $a\ge0$ ）、状態フィードバックを行って閉ループ系を漸近安定（ $a-bf<0$ ）にできることです。そのようなゲイン $f$ を決定する問題を状態フィードバックによる安定化問題と言います。

図1 状態フィードバックによる閉ループ系

[2] １次系(1)に対する安定化状態フィードバック(2)による閉ループ系に対して、評価関数

$\displaystyle{(4)\quad \boxed{J=\int_0^\infty(q^2x^2(t)+r^2u^2(t))\,dt} }$

を最小化するように $f$ を決めるLQ制御問題を考えます。

●閉ループ系における状態の振舞いは次式で与えられます。

$\displaystyle{(5)\quad x(t)=e^{(a-bf)t}x(0) }$

したがって、評価関数(4)の第１項は、次式のように計算されます。

$(6)\quad \begin{array}{lll} \displaystyle{J_x=\int_0^\infty q^2x^2(t)\,dt}\\ \displaystyle{=\int_0^\infty q^2e^{2(a-bf)t}x^2(0)\,dt}\\ \displaystyle{=q^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty}\\ \displaystyle{=\frac{q^2x^2(0)}{2(a-bf)}\left[\underbrace{e^{2(a-bf)\infty}}_{0}-\underbrace{e^{2(a-bf)0}}_{1}\right]}\\ \displaystyle{=-\frac{q^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0)} \end{array}$

これから $f$ を大きくすれば、 $J_x$ はいくらでも小さくできることが分かります。

●閉ループ系における入力の振舞い $u=-fx$ は次式で与えられます。

$\displaystyle{(7)\quad u(t)=-fe^{(a-bf)t}x(0) }$

したがって、評価関数(4)の第２項は、次式のように計算されます。

$(8)\quad \begin{array}{lll} \displaystyle{J_u=\int_0^\infty r^2u^2(t)\,dt}\\ \displaystyle{=\int_0^\infty r^2f^2e^{2(a-bf)t}x^2(0)\,dt}\\ \displaystyle{=r^2f^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty}\\ \displaystyle{=-\frac{r^2f^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0)} \end{array}$

$J_u$ が下に凸で、 $f=\frac{2b}{a}$ において最小値を取ることは、次式から分かります。

$(9)\quad \begin{array}{lll} \displaystyle{\frac{dJ_u}{df}=\frac{d}{df}\left(-\frac{r^2f^2}{2(a-bf)}x^2(0)\right)}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\left(\frac{2f}{(a-bf)}+(-1)\frac{f^2}{(a-bf)^2}(-b)\right)}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\frac{2(a-bf)f+f^2b}{(a-bf)^2}}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\frac{(2a-bf)f}{(a-bf)^2}} \end{array}$

$(10)\quad \begin{array}{lll} \displaystyle{\frac{d^2J_u}{df^2}=\frac{-r^2x^2(0)}{2}\left(\frac{2a-2bf}{(a-bf)^2}+(-2)\frac{(2a-bf)f}{2(a-bf)^3}(-b)\right)}\\ \displaystyle{=-r^2x^2(0)\left(\frac{1}{a-bf}+b\frac{2(a-bf)f+bf^2}{(a-bf)^3}\right)}\\ \displaystyle{=-r^2x^2(0)\frac{(a-bf)^2+2(a-bf)bf+b^2f^2}{(a-bf)^3}}r\\ \displaystyle{=-r^2x^2(0)\frac{(a-bf+bf)^2}{2(a-bf)^3}\}\\ \displaystyle{=-r^2x^2(0)\frac{a^2}{2(a-bf)^3}>0} \end{array}$

いま、 $b=1,q=1,r=1$ として、 $a=-1$ と $a=1$ の場合について、 $J_x$ 、 $J_u$ 、 $J=J_x+J_u$ の外形を描くと、次図のようになります。ここで、 $J$ には○で示す最小値が生まれていることに注意してください。

●それでは $J$ の最小値を与える $f$ を求めてみます。評価関数は

$(11)\quad \begin{array}{lll} \displaystyle{J=-\frac{q^2x^2(0)}{2(a-bf)}-\frac{r^2f^2x^2(0)}{2(a-bf)}}\\ \displaystyle{=\underbrace{-\frac{q^2+r^2f^2}{2(a-bf)}}_{\Pi}x^2(0)} \end{array} }$

のように書けますので、 $J$ の最小化は、次の $\Pi$ の最小化と等価です。

$\displaystyle{(12)\quad \boxed{\Pi=-\frac{q^2+r^2f^2}{2(a-bf)}} }$

図2 評価関数の概形

●この極値を求めるために $f$ で微分すると

$(13)\quad \begin{array}{lll} \displaystyle{\frac{d\Pi}{df}=-\frac{2r^2f}{2(a-bf)}-(-1)\frac{q^2+r^2f^2}{2(a-bf)^2}(-b)}\\ \displaystyle{=-\frac{2r^2(a-bf)f+(q^2+r^2f^2)b}{2(a-bf)^2}}\\ \displaystyle{=\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^2}} \end{array}$

これより

$\displaystyle{(14)\quad \frac{d\Pi}{df}=0 \ \Rightarrow \ bf^2-2af-\left(\frac{q}{r}\right)^2b=0 }$

ここで、 $f$ は $a-bf<0$ を満足させるものでなければならないので

$\displaystyle{(15)\quad \boxed{f=\frac{1}{b}\left(a+\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}\right)} }$

となります。そしてこれが一意に定まることは、以下のように $\Pi$ が下に凸でであることから分かります。

$(16)\quad \begin{array}{lll} \displaystyle{\frac{d^2\Pi}{df^2}=\frac{2br^2f-2ar^2}{2(a-bf)^2}+(-2)\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^3}(-b)}\\ \displaystyle{=\frac{-r^2}{(a-bf)}+b\frac{-2r^2(a-bf)f-br^2f^2-q^2b}{(a-bf)^3}}\\ \displaystyle{=\frac{-r^2(a-bf)^2-2r^2(a-bf)bf-r^2b^2f^2-q^2b^2}{(a-bf)^3}}\\ \displaystyle{=\frac{-r^2(a-bf+bf)^2-q^2b^2}{2(a-bf)^3}}\\ \displaystyle{=\frac{-r^2a^2-q^2b^2}{2(a-bf)^3}>0} \end{array}$

[2] LQ制御問題の解は求まったので、今度は評価関数(4)の重み係数 $q$ と $r$ の与え方について考えます。まず被積分項 $q^2x^2(t)+r^2u^2(t)$ は物理量の２乗和になっていますので、予め代表時間と代表長さによる無次元化を行っておく必要があります。たとえば、状態変数と入力変数がとりうる範囲を想定して

$\displaystyle{(17)\quad |x(t)|<M_x, |u(t)|<M_u\Rightarrow q=\frac{1}{M_x}, r=\frac{1}{M_u} }$

とすると、重み係数は物理単位の無次元化の役割を果たすことになります。次に(15)による閉ループ系を考えます。

$\displaystyle{(18)\quad \dot{x}(t)=(a-bf)x(t)=-\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}x(t) }$

ここで、 $a=\frac{1}{T},b=\frac{K}{T}$ の場合を考えると

$\displaystyle{(19)\quad \dot{x}(t)=-\frac{1}{T}\sqrt{1+\left(\frac{q}{r}\right)^2K^2}x(t) =-\frac{1}{\frac{T}{\sqrt{1+\left(\frac{q}{r}\right)^2K^2}}}x(t) }$

となって、新しい時定数は元の時定数 $T$ より小さくなることが分かります。これを決めるパラメータは重み係数の比 $q/r$ で、これが大きいほど時定数 $T$ は小さくなります。

●パラメータ $\rho$ を導入して重み係数を次のように設定します。

$\displaystyle{(20)\quad q=\frac{1}{M_x},\ r=\frac{1}{\rho M_u}\ \Rightarrow\ \frac{q}{r}=\rho\frac{M_u}{M_x}}$

ここで、 $\frac{M_u}{M_x}$ は固定されますので、設計パラメータは $\rho$ のみになります。

$a=-1$ 、 $b=1$ 、 $x(0)=1$ 、 $M_x=1$ 、 $M_u=1$ として、設計パラメータを $\rho=2,4,8$ と変えてみると、次の応答を得ます。

図3 重み係数の比 $\rho=q/r$ による応答の変化

演習A03-1…Flipped Classroom

$1^\circ$ 図2のグラフを描け。またパラメータ $\rho=q/r$ の値による状態と入力の振る舞いを観察し、トレードオフの関係について述べよ。

MATLAB

%a03_1.m
clear all, close all
a=-1; b=1; c=1;
Mx=1; Mu=1; q=1/Mx; r0=1/Mu
x0=1; 
t=0:0.1:4;
%-----
rho=2; r=r0/rho;
f1=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys1=ss(a-b*f1,b,[c;-f1],[]);
%-----
rho=4; r=r0/rho;
f2=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys2=ss(a-b*f2,b,[c;-f2],[]);
%-----
rho=8; r=r0/rho;
f3=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys3=ss(a-b*f3,b,[c;-f3],[]);
%-----
initial(sys1,sys2,sys3,x0,t)
grid
legend('rho=2','rho=4','rho=8'); 
title('LQ Control of 1st-order System')
xlabel('time')
ylabel('u(t) and x(t)')
%eof

SCILAB

coming soon

$2^\circ$ ハミルトン行列と呼ばれる

$\displaystyle{(21)\quad M=\left[\begin{array}{cc} a & -r^{-2}b^2 \\ -q^2 & -a \end{array}\right] }$

の安定な固有値

$\displaystyle{(22)\quad \lambda=-\sqrt{a^2+r^{-2}b^2q^2} }$

に対応する固有ベクトルは

$\displaystyle{(23)\quad \left[\begin{array}{cc} v_1 \\ v_2 \end{array}\right] =\left[\begin{array}{cc} 1 \\ \frac{-a-\sqrt{a^2+r^{-2}b^2q^2}}{-r^{-2}b^2} \end{array}\right] }$

となることを確かめ、次式を計算せよ。

$\displaystyle{(24)\quad f=r^{-2}bv_2v_1^{-1} }$

状態オブザーバ…Homework

[3] 状態フィードバックを使うには状態変数の変化をセンサで捉えてフィードバックする必要があるのですが、状態方程式が分かっているのですから、これを計算機内で解いてやれば良いのではないでしょうか？

いま制御対象（１次系）の状態空間表現を

$\displaystyle{(25)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t),\ x(0)\ne0\\ y(t)=cx(t) \end{array}\right. }$

とします。これを計算機で解きますが、その状態と出力を区別して

$\displaystyle{(26)\quad \left\{\begin{array}{lll} \dot{\hat{x}}(t)=a\hat{x}(t)+bu(t),\ \hat{x}(0)=0\\ \hat{y}(t)=c\hat{x}(t) \end{array}\right. }$

で表します。(26)から(25)の第１式どうしを辺々引くと

$\displaystyle{(27)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=a(\hat{x}(t)-x(t)),\ \hat{x}(0)\ne0 }$

を得ます。ここで $a<0$ でない限り

$\displaystyle{(28)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }$

とはなりません。そこで、図4に示すような工夫が行われました。

図4 状態オブザーバの考え方

● いま(2)に次のような補正項を加えます。

$\displaystyle{(29)\quad \dot{\hat{x}}(t)=a\hat{x}(t)+bu(t)-h(c\hat{x}(t)-y(t)),\ \hat{x}(0)=0 }$

すなわち

$\displaystyle{(30)\quad \dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+bu(t)+hy(t),\ \hat{x}(0)=0 }$

ここで、 $h$ は適当な定数です。(30)から(26)の第１式どうしを辺々引くと

$\displaystyle{(31)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=(a-hc)(\hat{x}(t)-x(t)),\ \hat{x}(0)\ne0 }$

を得ます。ここで $a-hc<0$ であるとすると

$\displaystyle{(32)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }$

が成り立ちます。すなわち状態オブザーバ(30)の状態 $\hat{x}(t)$ は、制御対象の状態 $x(t)$ を漸近的に推定することができます。このとき $a-hc<0$ を満たす(30)を状態オブザーバ、 $h$ をオブザーバゲインと呼びます。

● １次系

$\displaystyle{(33)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t),\ x(0)\ne0\\ y(t)=cx(t) \end{array}\right. }$

に対して、状態オブザーバ

$\displaystyle{(34)\quad \boxed{\dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+bu(t)+hy(t),\ \hat{x}(0)=0\quad (a-hc<0)} }$

によって推定された状態 $\hat{x}$ を用いる状態フィードバック

$\displaystyle{(35)\quad u(t)=-f\hat{x}(t) }$

による閉ループ系の状態方程式は、次式となります。

$\displaystyle{(36)\quad \left[\begin{array}{cc} \dot{x}(t)\\ \dot{\hat{x}}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a& -bf\\ hc& a-hc-bf \end{array}\right] }_{A_F} \left[\begin{array}{cc} x(t)\\ \hat{x}(t) \end{array}\right] }$

これから閉ループ系の安定性を判定したいのですが、このままでは見当がつきません。そこで、状態の推定誤差を $e(t)=\hat{x}(t)-x(t)$ とおくと、状態変数を変更した閉ループ系の状態方程式として、次式を得ます。

$\displaystyle{(37)\quad \boxed{ \left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a-bf & -bf\\ 0 & a-hc \end{array}\right] }_{A_F'} \left[\begin{array}{cc} x(t)\\ e(t) \end{array}\right]} }$

実際、状態の推定誤差の定義から、

$\displaystyle{(38)\quad \hat{x}(t)=x(t)+e(t) }$

に注意すると、状態フィードバック(35)は

$\displaystyle{(39)\quad u(t)=-fx(t)-fe(t) }$

これによる閉ループ系は、 $\dot{x}(t)=ax(t)+bu(t)$ に代入して

$\displaystyle{(40)\quad \dot{x}(t)=(a-bf)x(t)-bfe(t) }$

となります。これが(37)の第１式に他なりません。また、第２式は(31)そのもの、すなわち次式となります。

$\displaystyle{(41)\quad \dot{e}(t)=(a-hc)e(t) }$

閉ループ系の状態方程式は(40)と(41)に分離されており、 $a-bf$ と $a-hc$ は共に負となるように $f$ と $h$ を独立に設計してもよいことが分かります。この事実は分離定理と呼ばれています。

なお状態フィードバックを状態オブザーバの出力を用いて実施するコントローラ

$\displaystyle{(42)\quad \boxed{ \left\{\begin{array}{l} \dot{\hat{x}}(t)=(a-bf-hc)\hat{x}(t)+hy(t),\ \hat{x}(0)=0\\ u(t)=-f\hat{x}(t) \end{array}\right.} }$

は、オブザーバベースコントローラと呼ばれています。

演習A03-2…Flipped Classroom

$1^\circ$ １次系 $\dot{x}(t)=ax(t)+bu(t),y(t)=cx(t)$ に対する状態フィードバック $u(t)=-fx(t)$ と，状態オブザーバ $\dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+hy(t)+bu(t)$ （ $a-hc<0$ ）を考える。状態オブザーバによって推定された状態を用いる状態フィードバック $u(t)=-f\hat{x}(t)$ による閉ループ系の状態方程式を求めよ。

$2^\circ$ 閉ループ系の $A$ 行列の固有値が
$\lambda_1=-\frac{1}{T_1}$ ， $\lambda_2=-\frac{1}{T_2}$ （ $T_2>T_1>0$ ）となるように， $f$ と $h$ を定めよ。
$f$ と $h$ に， $\lambda_1$ と $\lambda_2$ をどう対応させたか，およびその理由を述べよ。

１次系のLQG制御…Homework

[5] １次系(33)に対して、オブザーバベースコントローラ(42)による閉ループ系に対して、評価関数

$\displaystyle{(43)\quad J=\int_0^\infty(q^2y\,^2(t)+r^2u\,^2(t))\,dt \label{eq7.2.7}}$

を最小化するように、状態フィードバックゲイン $f$ と状態オブザーバゲイン $h$ を決定するLQG制御問題を考えます。

●評価関数は次のように計算できます。

状態方程式

$\displaystyle{(44)\quad \underbrace{ \left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right] }_{\dot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{cc} a-bf & -bf\\ 0 & a-hc \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{cc} x(t) \\ e(t) \end{array}\right] }_{x_{CL}(t)} }$

と出力方程式

$\displaystyle{(45)\quad \underbrace{ \left[\begin{array}{cc} qx(t)\\ ru(t) \end{array}\right] }_{y_{CL}(t)} = \underbrace{ \left[\begin{array}{cc} q & 0 \\ -rf & -rf \end{array}\right] }_{C_{CL}} \underbrace{ \left[\begin{array}{cc} x(t) \\ e(t) \end{array}\right] }_{x_{CL}(t)} \label{eq7.2.5}}$

で表される閉ループ系の応答は

＜このパートは最初は読み飛ばして結構です＞
$\displaystyle{(46)\quad y_{CL}(t) = C_{CL}\exp(A_{CL}t) x_{CL}(0) \label{eq7.2.6}}$

このとき

$(47)\quad \begin{array}{lll} \displaystyle{J=\int_0^\infty y_{CL}^T(t)y_{CL}(t)\,dt}\\ \displaystyle{=\int_0^\infty x_{CL}^T(0)\exp(A_{CL}^Tt)C_{CL}^TC_{CL}\exp(A_{CL}t)x_{CL}(0)\,dt}\\ \displaystyle{=x_{CL}^T(0)\int_0^\infty \exp(A_{CL}^Tt)\underbrace{C_{CL}^TC_{CL}}_{Q_{CL}}\exp(A_{CL}t)\,dt\,x_{CL}(0)}\\ \displaystyle{=x_{CL}^T(0)\underbrace{\int_0^\infty \exp(A_{CL}^Tt)Q_{CL}\exp(A_{CL}t)\,dt}_{\Pi}\,x_{CL}(0)}\\ \displaystyle{=x_{CL}^T(0)\Pi x_{CL}(0)} \end{array}$

ここで、初期値

$\displaystyle{(48)\quad x_{CL}(0)= \left[\begin{array}{cc} x(0)\\ e(0) \end{array}\right]= \left[\begin{array}{cc} x(0)\\ \hat{x}(0)-x(0) \end{array}\right] \label{eq7.2.9}}$

をどう選ぶかですが、図4に示す閉ループ系を考えます。すなわち、(33)の代わりに、

$\displaystyle{(49)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t)+\sigma n_x(t),\ x(0)\ne0\\ y(t)=cx(t)+\rho n_y(t) \end{array}\right. \label{eq7.2.10}}$

を考えます。ここで、新たに加えられた入力 $n_x(t)$ と $n_y(t)$ はそれぞれシステムノイズと観測ノイズとよばれ、定常正規確率過程として扱うのが普通です。特に、 $\sigma^2$ と $\rho^2$ は分散を表しています。ここでは確率的な議論を避けるため、 $n_x(t)$ と $n_y(t)$ は確定的なインパルス入力であると仮定します。

図5 オブザーバベース・コントローラによる閉ループ系の評価

このとき、閉ループ系の状態方程式は、次式となります。

$\displaystyle{(50)\quad &&\left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a-bf & -bf \\ 0 & a-hc \end{array}\right] }_{A_{CL}} \left[\begin{array}{cc} x(t)\\ e(t) \end{array}\right]\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} \sigma & 0\\ -\sigma & h\rho \end{array}\right] }_{B_{CL}} \left[\begin{array}{cc} n_x(t)\\ n_y(t) \end{array}\right] \label{eq7.2.11}}$

出力方程式(45)の場合のインパルス応答は、インパルス入力を $n_x$ に与える場合と $x_y$ に与える場合の２通り

＜このパートは最初は読み飛ばして結構です＞
$\displaystyle{(51)\quad y^{(1)}_{CL}(t) = C_{CL}\exp(A_{CL}t) \underbrace{ \left[\begin{array}{cc} \sigma \\ -\sigma \end{array}\right] }_{x_{CL}(0)=B^{(1)}_{CL}} \label{eq7.2.12}}$

$\displaystyle{(52)\quad y^{(2)}_{CL}(t) = C_{CL}\exp(A_{CL}t) \underbrace{ \left[\begin{array}{cc} 0 \\ h\rho \end{array}\right] }_{x_{CL}(0)=B^{(2)}_{CL}} \label{eq7.2.13}}$

があるので、評価関数を改めて次式のように両者の２乗和とします。

$\displaystyle{(53)\quad J'=B_{CL}^{(1)T}\Pi B_{CL}^{(1)}+B_{CL}^{(2)T}\Pi B_{CL}^{(2)} \label{eq7.2.14}}$

これは、 $\Pi= \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]$ とおくと

$\displaystyle{(54)\quad J'=\pi_1\sigma^2+\pi_2(\sigma^2+h^2\rho^2)-2\pi_3\sigma^2 \label{eq7.2.15}}$

と計算できます。ここで、 $\pi_1,\pi_2,\pi_3$ を直接計算することは容易ではないのですが、

＜このパートは最初は読み飛ばして結構です＞
$A_{CL}$ が安定行列であることは， $\pi_1,\pi_2,\pi_3$ が，リャプノフ方程式

$\displaystyle{(55)\quad \begin{array}{lll} && \underbrace{ \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] }_{\Pi} \underbrace{ \left[\begin{array}{cc} a-bf & -bf \\ 0 & a-hc \end{array}\right] }_{A_{CL}}\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} a-bf & 0 \\ -bf & a-hc \end{array}\right] }_{A_{CL}^T} \underbrace{ \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] }_{\Pi}\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} q^2+r^2f^2 & r^2f^2 \\ r^2f^2 & r^2f^2 \end{array}\right] }_{C_{CL}^TC_{CL}>0} = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array} \label{eq7.2.16}}$

すなわち

$\displaystyle{(56)\quad 2(a-bf)\pi_1+q^2+r^2f^2=0 \label{eq7.2.17}}$
$\displaystyle{(57)\quad -bf\pi_1+(2a-bf-hc)\pi_3+r^2f^2=0 \label{eq7.2.18}}$
$\displaystyle{(58)\quad 2(a-hc)\pi_2-2bf\pi_3+r^2f^2=0 \label{eq7.2.19}}$

の解であることに等しいという結果を用います。

(56)より $\pi_1$ を求めれば，評価関数は

$\displaystyle{(59)\quad J'= \underbrace{ \frac{q^2+r^2f^2}{-2(a-bf)} }_{\pi_1} \sigma^2 +\pi_2(\sigma^2+h^2\rho^2)-2\pi_3\sigma^2 \label{eq7.2.20}}$

となります。ここで，第1項は，(11)で， $x(0)=\sigma$ としたものであることに注意します。いま， $\hat{x}(0)=x(0)$ かつ $w(t)=v(t)=0$ という場合を考えてみると，(50)の第2行目より

$\displaystyle{(60)\quad \dot{e}(t)=(a-hc)e(t),e(0)=0\ \Rightarrow\ e(t)=0\ \Rightarrow\ \hat{x}(t)=x(t) \label{eq7.2.21}}$

だから，(35)は状態フィードバック $u(t)=-fx(t)$ に相当します。この場合の議論はすでに行っており， $\pi_2,\pi_3$ は存在しないので，評価関数(59)は第1項だけとなり，これを最小化する $f$ は，リッカチ方程式

$\displaystyle{(61)\quad -r^{-2}b^2\pi_1^2+2a\pi_1+q^2=0 \label{eq7.2.22}}$

の解 $\pi_1>0$ を用いて次のように与えられます。

$\displaystyle{(62)\quad f=r^{-2}b\pi_1 \label{eq7.2.23}}$

そうすると，評価関数(59)の第2項と第3項は，状態オブザーバを用いるために生じると考えられます。このように，理想的な状態フィードバックに比べて評価関数の値が増えることを評価関数の劣化と呼びます。

つぎに，仮定 $f=r^{-2}b\pi_1$ ， $a-bf<0$ ， $a-hc<0$ のもとでは，(57)より $\pi_3=0$ を得ます。このとき，式(58)は

$\displaystyle{(63)\quad 2(a-hc)\pi_2+r^2f^2=0 \label{eq7.2.24}}$

となり，これから $\pi_2$ が求まります。したがって，評価関数(59)は

$\displaystyle{(64)\quad J'= \underbrace{ \frac{q^2+r^2f^2}{-2(a-bf)} }_{\#1} \sigma^2 + \underbrace{ \frac{\sigma^2+h^2\rho^2}{-2(a-hc)} }_{\#2} r^2f^2 \label{eq7.2.25}}$

と書けます。評価関数の劣化分を最小化するためには，#2を $h$ について最小化すればよいことになります。幸いなことに，#2は#1と同形であるので，これを最小化する $h$ の決定法は $f$ と同様に行ないます。

●以上の議論をまとめると、次のようになります。

アルゴリズム＜1次系のLQG制御＞

入力パラメータ
$a,\,b,\,c,\,q>0,\,r>0,\,\sigma>0,\,\rho>0$
出力パラメータ
$a_K,\,b_K,\,c_K$

ステップ1: 状態フィードバックゲイン $f$ の決定

$\displaystyle{(65)\quad \boxed{{\rm\bf CARE:\ }\displaystyle{-\frac{1}{r^2}b^2\Pi^2+2a\Pi+q^2=0}} }$

を解いて， $\Pi>0$ を求め， $f$ をつぎのように定める。

$\displaystyle{(66)\quad \boxed{f=\frac{1}{r^2}b\Pi} }$

ステップ2: オブザーバゲイン $h$ の決定

$\displaystyle{(67)\quad \boxed{{\rm\bf FARE:\ }\displaystyle{-\frac{1}{\rho^2}c^2\Gamma^2+2a\Gamma+\sigma^2=0}} }$

を解いて， $\Gamma>0$ を求め， $h$ をつぎのように定める。

$\displaystyle{(68)\quad \boxed{h=\frac{1}{\rho^2}c\Gamma} }$

ステップ3: オブザーバベースコントローラの構成

$\displaystyle{(69)\quad \left\{\begin{array}{l} \dot{x}_K(t)=a_Kx_K(t)+b_Ky(t)\\ u(t)=c_Kx_k(t) \end{array}\right. }$

ただし

$\displaystyle{(70)\quad a_K=a-hc-bf,\ b_K=h,\ c_K=-f }$

● $a=-1$ 、 $b=1$ 、 $c=1$ 、 $q=1$ 、 $r=1$ 、 $\sigma=1$ 、 $\rho=0.2$ 、 $x(0)=1$ 、 $\hat{x}(0)=0$ として、状態フィードバックだけとオブザーバベースコントローラを用いた場合とを比較してみると、次の応答を得ます。

図6 状態フィードバックとオブザーバベースコントローラを用いた場合の応答の変化

演習A03-3…Flipped Classroom

$1^\circ$ 図6のグラフを描け。

MATLAB

%a03_2.m
 clear all, close all
 a=-1; b=1; c=1; 
 q=1; r=1; 
 sig=1; rho=0.2;
 f=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
 sys1=ss(a-b*f,[],[c;-f],[]);
 x0=1; 
 t=0:0.1:4;
 initial(sys1,x0,t), hold on
 h=1/c*(a+sqrt(a^2+(sig/rho)^2*c^2));
 AF=[a -b*f;
     h*c a-h*c-b*f];
 sys2=ss(AF,[],diag([c,-f]),[]);
 initial(sys2,[x0;0],t)  
 grid
 legend('sfc','obc')
 title('LQ OBC of 1st-order System')
 xlabel('time')
 ylabel('u(t) and x(t)')
%eof

SCILAB

coming soon

$2^\circ$ 上のシミュレーションにおいて、 $a-hc<a-bc$ となっていることを確かめよ。

Note A03　ラグランジュの未定定数法

１入力１出力１次系に対する最適制御（LQ制御）を考えましたが、一般の多入力多出力高次系の場合も同様な議論ができて、大変有用な制御手段として知られています。その説明はあとで行ないますが、そこではラグランジュの未定定数法によるアプローチが採用されます。

(11)の評価関数 $J$ を最小化する問題を、 $\Gamma$ をラグランジュの未定定数として

$\displaystyle{(1)\quad J=\Pi x^2(0)+\Gamma(2(a-bf)\Pi+q^2+r^2f^2) }$

を最小化する問題として考え直してみます。ここで、制約条件は、リャプノフ方程式と呼ばれる

$\displaystyle{(2)\quad \boxed{2(a-bf)\Pi+q^2+r^2f^2=0} }$

ですが、 $\Pi>0$ は $a-bf<0$ を意味することに注意します。そこで、必要条件として次を得ます。

$\displaystyle{(3)\quad \frac{\partial J}{\partial \Pi}&=&x^2(0)+2(a-bf)\Gamma=0\Rightarrow\Gamma>0 }$
$\displaystyle{(4)\quad \frac{\partial J}{\partial f}&=&(-2b\Pi+2r^2f)\Gamma=0\Rightarrow f=r^{-2}b\Pi }$
$\displaystyle{(5)\quad \frac{\partial J}{\partial \Gamma}&=&2(a-bf)\Pi+q^2+r^2f^2=0 }$

ここで、第２式から得られる $f=r^{-2}b\Pi$ を第３式に代入して

$\displaystyle{(6)\quad 2(a-br^{-2}b\Pi)\Pi+q^2+r^2r^{-4}b^2\Pi^2=0 }$

すなわち、リッカチ方程式と呼ばれる $\Pi$ の２次方程式

$\displaystyle{(7)\quad 2a\Pi-r^{-2}b^2\Pi^2+q^2=0 }$

を得ます。これから $\Pi>0$ を

$\displaystyle{(8)\quad \Pi=\frac{a+\sqrt{a^2+r^{-2}q^2b^2}}{r^{-2}b^2} }$

を求めて、 $F$ は

$\displaystyle{(9)\quad f=r^{-2}b\Pi=\frac{1}{b}\left(a+\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}\right) }$

のように得られます。十分性の議論はここでは省きます。

１次系の時間応答

投稿日時: 2021年10月8日投稿者: cacsd

１次系の時間応答…Homework

[1] 次の１入力１出力１次系の状態空間表現を考えます。

$\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}(t)=ax(t)+bu(t)\\ y(t)=cx(t) \end{array}\right. }$

これから、１次系の時間応答の表現式

$\displaystyle{(2)\quad y(t)=ce^{at}x(0)+\int_0^tce^{a(t-\tau)}bu(\tau)d\tau}$

を得ます（Note A02-1参照）。ここで、インパルス応答と呼ばれる

$\displaystyle{(3)\quad \boxed{g(t)=ce^{at}b}}$

を定義しますと、上式は

$\displaystyle{(4)\quad \boxed{y(t)=\underbrace{ce^{at}x(0)}_{zero-input\ response} +\underbrace{\int_0^tg(t-\tau)u(\tau)d\tau}_{zero-state\ response}} }$

すなわち、１次系の時間応答は、零入力応答（第１項）と零状態応答（第２項）の和でとなります。

●実際、 $a=-1,b=1,c=1$ 、 $x(0)=1$ 、 $u(t)=\sin(10t)$ について $y(t)$ を求めてみると次のグラフが得られます。

図1 どれが零状態応答か？

[2] (1)への入力 $u$ として次図の時間関数を考えます。

図2 ステップ入力とインパルス入力( $T_s\rightarrow0$ )

左図において $T_s\rightarrow0$ としたものがステップ入力、

$\displaystyle{(5)\quad u(t)= \left\{\begin{array}{ll} 0 & (t\le0)\\ 1 & (t>0) \end{array}\right.}$

です。これを(2)に代入すると、

$\displaystyle{(6)\quad y(t)=\int_0^tg(t-\tau)\cdot 1 d\tau=\int_{t}^0g(t')(-dt')=\int_0^tg(t')dt'=\int_0^tg(\tau)d\tau}$

より、ステップ応答が次のように表されます。

$\displaystyle{(7)\quad \boxed{y(t)=\int_0^tg(\tau)d\tau}}$

すなわち、インパルス応答を積分してステップ応答が得られます。逆に、ステップ応答を微分してインパルス応答が得られることも分かります。

図2の右図は左図の微分 $\dot{u}$ を表していますが、 $T_s\rightarrow0$ としたものがインパルス入力（デルタ関数）です。実は(3)で定義したインパルス応答はこのインパルス入力に対する応答です。このことを理解するために、いま１次系の状態方程式を微分して

$\displaystyle{(8)\quad \ddot{x}(t)=a\dot{x}(t)+b\dot{u}(t)}$

を得ます。これは１次系に入力 $\dot{u}$ を与えた場合の応答は $\dot{x}(t)$ であることを示しています。図2の左図で $T_s\rightarrow0$ としたステップ入力に対する応答がステップ応答です。図2の右図で $T_s\rightarrow0$ としたインパルス入力はステップ入力の微分ですから、その応答はステップ応答の微分として得られ、インパルス応答となります。すなわち、ステップ応答を微分するとインパルス応答が得られることが分かります。

[3] 特別な３種類の入力（ステップ入力、インパルス入力、正弦波入力）に対する零状態応答（ステップ応答、インパルス応答、正弦波応答）を考え、これらの相互関係を調べます。

図3 １次系の時間応答の相互関係

一般に、「システムを知る」ということは「どんな入力を与えられても出力を推測できる」ことを意味します。したがって、平衡状態（ $x(0)=0$ ）にある場合、インパルス応答さえ分かれば、畳み込み積分を行って出力（零状態応答）が計算できます。この意味で、インパルス応答は重要なのです。また、インパルス応答自体の計算は、 $x(0)=b$ としたときの零入力応答を求めて行います。ステップ応答が分かれば、これを微分してインパルス応答が分かります。また、インパルス応答をラプラス変換したものが、正弦波応答と密接に関係しています。このことはあとで詳しく述べます。

[4] 次のようにパラメータライズされた漸近安定な１次系を考えます。 $T$ 、 $K$ はそれぞれ時定数、定常ゲインと呼ばれます。

$\displaystyle{(9)\quad \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t)\\ y(t)=x(t) \end{array}\right. }$

このインパルス応答は次式で与えられます。

$\displaystyle{(10)\quad g(t)=\frac{K}{T}e^{-\frac{1}{T}t}}$

これを積分して

$\displaystyle{(11)\quad y(t)=\frac{K}{T}\int_0^te^{-\frac{1}{T}\tau}d\tau=\frac{K}{T}\left[\frac{e^{-\frac{1}{T}\tau}}{-\frac{1}{T}}\right]_0^t =-K(e^{-\frac{1}{T}t}-1)=K(1-e^{-\frac{1}{T}t}) }$

すなわちステップ応答は次式で表されます。

$\displaystyle{(12)\quad \boxed{y(t)=K(1-e^{-\frac{1}{T}t})}}$

ここで、 $t=T$ のときは

$\displaystyle{(13)\quad y(T)=K(1-\frac{1}{e})=0.632K}$

となります。また、 $t=0$ での接線は

$\displaystyle{(14)\quad y(t)=\dot{x}(0)t=\frac{K}{T}t}$

ですから、 $x(T)=K$ となります。この時定数 $T$ は、最終値 $K$ の63.2%まで到達する時間、または初期速度で到達した場合の時間（初期時刻における接線が $K$ と交わる時刻）として特徴付けられます。

●実際、 $T=0.5,1,2$ 、 $K=1$ について、ステップ応答 $x(t)$ を求めてみると次のグラフが得られます。最終値 $K=1$ に到達する時間は大体 $5T$ であることが分かります。

図4 １次系のステップ応答から時定数を求める補助線とは？

[5] 時定数は平衡状態に復帰する時間の目安にもなります。漸近安定な１次自由系を次のようにパラメータライズします。

$\displaystyle{(15)\quad \dot{x}(t)=-\frac{1}{T}x(t)}$

この解は

$\displaystyle{(16)\quad x(t)=e^{-\frac{1}{T}t}x(0)}$

ですから、 $t=T$ のときは

$\displaystyle{(17)\quad x(T)=\frac{1}{e}x(0)=0.367x(0)}$

となります。また、 $t=0$ での接線は

$\displaystyle{(18)\quad x(t)=\dot{x}(0)t+x(0)=-\frac{1}{T}x(0)t+x(0)}$

ですから、 $x(T)=0$ となります。この時定数 $T$ は、初期値の36.7%まで復帰する時間、または初期速度で復帰した場合の時間（初期時刻における接線が0と交わる時刻）として特徴付けられます。

●実際、時刻 $t=0$ において平衡状態 $x=0$ が乱され $x(0)=1$ になったとし、 $T=0.5,1,2$ について $x(t)$ を求めてみると次のグラフが得られます。平衡状態に復帰する時間は大体 $5T$ であることを表します。

<
図5 １次自由系において平衡状態に復帰する時間（時定数）を定める補助線とは？

演習A02-1…Flipped Classroom

$1^\circ$ 図1の時間応答を描くプログラムを作成せよ。

MATLAB


%a01_1.m
clear all, close all
a=-1; b=1; c=1;
sys=ss(a,b,c,[]);
t=0:0.01:10;
u=sin(10*t);
lsim(sys,u,t,1), hold on
lsim(sys,u,t,0)
lsim(sys,u*0,t,1)
grid
title('Time Response of 1st-order System')
xlabel('time')
ylabel('x(t)')
legend('total resp','zero-state resp','zero-input resp')
%eof

SCILAB

coming soon

$2^\circ$ 図4のステップ応答#1を描き、適当な補助線を引いて、交点の座標から時定数を求めるプログラムを作成せよ。

MATLAB

%a02_2.m
clear all, close all
T1=0.5; T2=1; T3=2; K=1;
a1=-1/T1; a2=-1/T2; a3=-1/T3;
b1=K/T1; b2=K/T2; b3=K/T3; 
c=1;
sys1=ss(a1,b1,c,[]);
sys2=ss(a2,b2,c,[]);
sys3=ss(a3,b3,c,[]);
t=0:0.01:5;
step(sys1,sys2,sys3,t)
grid
hold on
plot([0,5],(1-1/exp(1))*[1 1])
T=ginput(3)
title('Step Responses of 1st-order System')
xlabel('time')
ylabel('x(t)')
legend('T1=0.5','T2=1','T3=2','1-1/e')
%eof

SCILAB

coming soon

[6] １次系(1)をラプラス変換（Note A02-3参照）して、 $x(0)=0$ とすると（零状態応答を考えているので）

$\displaystyle{(19)\quad s\hat{x}(s)-\underbrace{x(0)}_0=a\hat{x}(s)+b\hat{u}(s) \Rightarrow \hat{x}(s)=(s-a)^{-1}b\hat{u}(s) }$

$\displaystyle{(20)\quad \hat{y}(s)=c\hat{x}(s) \Rightarrow \hat{y}(s)=\underbrace{c(s-a)^{-1}b}_{\hat{G}(s)}\hat{u}(s) }$

ここで現れた

$\displaystyle{(21)\quad \boxed{\hat{g}(s)=c(s-a)^{-1}b}}$

を伝達関数と呼びます。伝達関数 $G(s)$ に $s=j\omega$ を代入した $G(j\omega)$ を周波数伝達関数と呼び、正弦波入力に対する時間応答、すなわち周波数応答を特徴づけることができます。

●漸近安定な１次系(9)の伝達関数は次式で与えられます。

$\displaystyle{(22)\quad \hat{g}(s)=\frac{K}{Ts+1} }$

演習A02-2...Flipped Classroom

$1^\circ$ １次系のステップ応答を次式を逆ラプラス変換して求めなさい。

$\displaystyle{(23)\quad \hat{y}(s)=\underbrace{\frac{K}{Ts+1}}_{\hat{g}(s)}\frac{1}{s}=K(\frac{1}{s}-\frac{T}{Ts+1}) }$

$2^\circ$ 次のような無駄時間と呼ばれる特殊な入出力関係を考えます。

$\displaystyle{(24)\quad y(t)=u(t-L)\quad(L>0) }$

これをラプラス変換して無駄時間の伝達関数は次式をとなることを示しなさい。

$\displaystyle{(25)\quad \boxed{\hat{g}(s)=e^{-sL}} }$

ヒント：

$(26)\quad \begin{array}{l} \displaystyle{\hat{y}(s)=\int_0^\infty u(t-L)e^{-st}\,dt =\int_{-L}^\infty u(\tau)e^{-s(\tau+L)}\,d\tau}\\ \displaystyle{=\underbrace{e^{-sL}}_{\hat{g}(s)}\underbrace{\int_0^\infty u(\tau)e^{-s\tau}\,d\tau}_{\hat{u}(s)}} \end{array}$

[7] 漸近安定な１次系(9)を考えます。正弦波入力に対する零状態応答すなわち周波数応答は、Note A02-2を参照して

$\displaystyle{(27)\quad y(t)=\frac{K}{\sqrt{1+(\omega T)^2}}(\sin(\omega t-\tan^{-1}\omega T)+e^{-\frac{1}{T}t}\sin(\tan^{-1}\omega T))}$

ここで、 $t\rightarrow\infty$ とすると

$\displaystyle{(28)\quad \boxed{y(t)\simeq\underbrace{\frac{K}{\sqrt{1+(\omega T)^2}}}_{|\hat{g}(j\omega)|} \sin(\omega t\underbrace{-\tan^{-1}\omega T}_{\angle\hat{g}(j\omega)})} }$
ただし
$\displaystyle{(29)\quad \hat{g}(s)=\frac{K}{Ts+1} }$

これは、入力が正弦波のときは、時間が十分立てば、出力も正弦波となることを示しています。その振幅と位相はそれぞれ $\hat{g}(j\omega)$ の絶対値と偏角となっています。

$|\hat{g}(j\omega)|$ をゲイン、 $\angle\hat{g}(j\omega)$ を位相と呼びます。このゲイン線図と位相線図をペアにして片対数グラフにプロットしたものをボード線図と呼びます。ゲインはdb値( $20{\log_{10}|\hat{g}(j\omega)|$ )、位相はdeg値( $\frac{180}{\pi}\angle\hat{g}(j\omega)$ )です。

●実際、 $T=0.5,1,2$ について、ボード線図を描いてみると次のグラフが得られます

図6 １次系のボード線図から時定数を求める補助線とは？

ちなみに

$\displaystyle{(30)\quad |\hat{G}(j\omega_b)|=\frac{|\hat{G}(j0)|}{\sqrt{2}}}$

を満足する $\omega_b$ は

$\displaystyle{(31)\quad \omega_b=\frac{1}{T}}$

となり、帯域幅と呼びます。実際、上図で-3dBの水平線との交点が $1/T$ であることが確かめられます。

演習A02-3...Flipped Classroom

$1^\circ$ 図1のゲイン曲線#1を描き、3dB下がる点の座をマウスをクリックして取得し、時定数を求めるプログラムを作成せよ。

MATLAB

%a02_3.m
clear all, close all
T1=0.5; T2=1; T3=2; K=1;
a1=-1/T1; a2=-1/T2; a3=-1/T3;
b1=K/T1; b2=K/T2; b3=K/T3; 
c=1;
sys1=ss(a1,b1,c,[]);
sys2=ss(a2,b2,c,[]);
sys3=ss(a3,b3,c,[]);
sys0=ss([],[],[],1/sqrt(2));
w=logspace(-1,1);
bode(sys1,sys2,sys3,sys0,w)
grid
W=ginput(3)
title('Bode Diagrams of 1st-order System')
xlabel('Freq')
ylabel('Gain')
%ylabel('Phase')
legend('T1=0.5','T2=1','T3=2','-3dB')
%eof

SCILAB

coming soon

Note A02-1　状態方程式の解

●(１次自由系を表す)微分方程式

$\displaystyle{(1)\quad \dot{x}(t)=ax(t)}$

の解は

$\displaystyle{(2)\quad x(t)=e^{at}x(0)}$

と表されます。これが解であることは元の微分方程式に代入すればすぐに確かめられ、また次のようにして導くことができます。

元の微分方程式を

$\displaystyle{(3)\quad \dot{x}(t)-ax(t)=0}$

と書いて、左から積分因数と呼ばれる $e^{-at}$ をかけると

$\displaystyle{(4)\quad e^{-at}\dot{x}(t)-ae^{-at}x(t)=0}$

すなわち

$\displaystyle{(5)\quad \frac{d}{dt}(e^{-at}x(t))=0}$

したがって、 $c$ を定数として

$\displaystyle{(6)\quad e^{-at}x(t)=c\ \Rightarrow\ x(t)=e^{at}c}$

ここで、 $t=0$ とおくと $c$ は初期値 $x(0)$ に等しいので

$\displaystyle{(7)\quad \boxed{x(t)=e^{at}x(0)}}$

と表されます。

●(１次系の状態方程式を表す)微分方程式

$\displaystyle{(8)\quad \dot{x}(t)=ax(t)+bu(t)}$

の解を求めてみます。上と同様にして

$\displaystyle{(9)\quad \frac{d}{dt}(e^{-at}x(t))=e^{-at}bu(t)}$

までは問題ないと思います。これを $0$ から $t'$ まで積分して

$\displaystyle{(10)\quad \left[e^{-at}x(t)\right]_0^{t'}=\int_0^{t'}e^{-at}bu(t)dt}$

すなわち

$\displaystyle{(11)\quad e^{-at'}x(t')-x(0)=\int_0^{t'}e^{-at}bu(t)dt}$

したがって

$\displaystyle{(12)\quad x(t')=e^{at'}x(0)+e^{at'}\int_0^{t'}e^{-at}bu(t)dt}$

ここで、 $t$ を $\tau$ に、 $t'$ を $t$ と置き換えれば

$\displaystyle{(13)\quad \boxed{x(t)=e^{at}x(0)+\int_0^te^{a(t-\tau)}bu(\tau)d\tau}}$

が得られます。

Note A02-2　正弦波応答

(１次系の状態方程式を表す)微分方程式

$\displaystyle{(1)\quad \dot{x}(t)=ax(t)+bu(t) }$

に対して、正弦波入力

$\displaystyle{(2)\quad u(t)=\sin\omega t}$

に対する解を求めます。これは公式

$\displaystyle{(3)\quad \int e^{ax}\sin bx\,dx=\frac{e^{ax}}{a^2+b^2}(a\sin bx-b \cos bx)}$

を用いて、次のように得られます。

$(4)\quad \begin{array}{l} \displaystyle{x(t)=\int_0^te^{a(t-\tau)}b \sin\omega\tau\,d\tau=e^{at}b\int_0^te^{-a\tau} \sin\omega\tau\,d\tau}\\ \displaystyle{=e^{at}b\left[\frac{e^{-a\tau}}{a^2+\omega^2}(-a\sin\omega\tau-\omega\cos\omega\tau)\right]_0^t}\\ \displaystyle{=\frac{e^{at}b}{a^2+\omega^2}(e^{-at}(-a\sin\omega t-\omega\cos\omega t)+\omega)}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin\omega t\frac{-a}{\sqrt{a^2+\omega^2}}-\cos\omega t\frac{\omega}{\sqrt{a^2+\omega^2}}+e^{+at}\frac{\omega}{\sqrt{a^2+\omega^2}})}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin\omega t\cos\theta-\cos\omega t\sin\theta+e^{+at}\sin\theta)}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin(\omega t-\theta)+e^{+at}\sin\theta)\quad (\theta=\tan^{-1}\frac{\omega}{-a})}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin(\omega t+\phi)-e^{+at}\sin\phi)\quad (\phi=-\theta)} \end{array}$

Note A02-3　ラプラス変換

時間関数 $x(t)$ のラプラス変換は次式によって定義されます。

$\displaystyle{(1)\quad \boxed{\hat{x}(s)=\int_0^\infty x(t)e^{-st}\,dt}}$

Note A02-4　dB値 $20{\log_{10}(\cdot)$

１次系の状態空間表現

投稿日時: 2021年10月8日投稿者: cacsd

一般には制御対象は非線形系ですが、ここでは重ね合わせの原理の成り立つ線形系（linear system）を考えます。また１入力１出力１次系（SISO 1st-order system）に限って、制御理論の概要を漸近安定性、時間応答、安定化制御、追値制御の４つに分けて述べます。

１次系の状態空間表現…Homework

[1] 　制御対象が１入力１出力１次系の場合、そのモデルとして、次の状態空間表現が用いられます。

$\displaystyle{(1)\quad \boxed{\left\{\begin{array}{ll} \dot{x}(t)=ax(t)+bu(t)&(x(t),u(t)\in{\bf R})\\ y(t)=cx(t)&(y(t)\in{\bf R}) \end{array}\right.} }$

第１式は状態方程式とよばれ、 $x(t),u(t)$ はそれぞれ時刻 $t$ における状態変数、入力変数（アクチュエータを操作する）を表しています。また、第２式は出力方程式とよばれ、 $y(t)$ は時刻 $t$ における出力変数（センサを用いて観測される）を表しています。 ${\bf R}$ は実数の集合を表します。以下では状態空間表現(1)を簡単に１次系と参照します。

この状態空間表現は図3のようなブロック線図で表されます。

図1 １入力１出力１次系状態空間表現のブロック線図

●具体例を示すために、次図のような走行する車を考えます。

図2 走行する車

ニュートンの運動第２法則「質量×加速度＝外力」を適用するため、車は質点とみなし、直線運動をしているとします。車の質量を $m$ 、時刻 $t$ における車の速度を $v(t)$ 、車に働く外力（制動力）を $f(t)$ とすると、車の運動方程式は次の微分方程式で与えられます。

$\displaystyle{(2)\quad m\dot{v}(t)=f(t)\qquad(v(0)\ne0) }$

ここで、 $\dot{v}(t)=\frac{d}{dt}v(t)$ は車の加速度を表します。もし空気による抗力 $cv(t)$ を考慮する場合は、車の運動方程式は次式となります。

$\displaystyle{(3)\quad m\dot{v}(t)=-cv(t)+f(t) }$

いま車は一定速度 $v^*$ で等速運動をしているとします。このとき加速度は零となるので

$\displaystyle{(4)\quad 0=-cv^*+f^* }$

を満たす一定の制動力 $f^*=cv^*$ が定まります。実際には、この等速運動が何らかの原因により乱されたとき、これを速やかに元に戻すことが要求されます。そこで(3)から(4)を辺々引き算すると次式を得ます。

$\displaystyle{(5)\quad \underbrace{\frac{d}{dt}(v(t)-v^*)}_{\dot{x}(t)}=\underbrace{-\frac{c}{m}}_{a}\underbrace{(v(t)-v^*)}_{x(t)}+\underbrace{\frac{1}{m}}_{b}\underbrace{(f(t)-f^*)}_{u(t)} }$

これが入力１出力１次系の状態方程式の一例です。ここで、 $x=0$ と $u=0$ が、平衡状態 $v^*$ とこれを維持する平衡入力 $f^*$ を表していることに留意します。

１次系の漸近安定性…Homework

[2]　1次系(1)で表される制御対象は平衡状態にあるとします（ $x=0, u=0$ ）。いま平衡状態が乱され、時刻 $t=0$ において $x(0)\ne 0$ となったとします。特に制動力を変えないとすると $u=0$ です。このとき、元の平衡状態 $x=0$ に復帰できるならば、１次系(1)は漸近安定(asymptotically stable)であると言います。そこで

$\displaystyle{(6)\quad \boxed{\dot{x}(t)=ax(t)\qquad(x(0)\ne0)} }$

の解

$\displaystyle{(7)\quad x(t)=e^{at}x(0) }$

の振舞いを調べてみます。 $t\rightarrow\infty$ のとき

$\displaystyle{(8)\quad x(t) \rightarrow \left\{\begin{array}{ll} 0&(a<0)\\ x(0)&(a=0)\\ \infty&(a>0) \end{array}\right. }$

となるので、１次系(1)の漸近安定性の必要十分条件は

$\displaystyle{(9)\quad \boxed{a<0} }$

であることがわかります。

このように漸近安定性は、(1)において $u=0$ とした自由系(unforced system)(6)に基づいて判定されることに留意してください。

●上の等速運動をしている車の例を考えます。いま突風による急激な速度変化のため $v(0)\ne v^*$ となり、特に制動力を変えないとすると $u=0$ なので

$\displaystyle{(10)\quad \underbrace{\frac{d}{dt}(v(t)-v^*)}_{\dot{x}(t)}=\underbrace{-\frac{c}{m}}_{a}\underbrace{(v(t)-v^*)}_{x(t)} \qquad(x(0)=v(0)-v^*\ne0) }$

を解いて、元の速度 $v^*$ に復帰できるかを判断できます。ここで、 $a<0$ ですから、漸近安定であることが確かめられます。

以上の議論は数学的には明らかですが、物理的にはどのような力が働いて平衡状態に戻るのでしょうか？それは空気抵抗 $-cv(t)$ が速度の増減によって、逆方向に働くからといえます。速度に比例する抗力を減衰力、位置に比例する抗力を復元力といいます。モノつくりにおいて、これらを内在させることをパッシブ制御、コントローラ（補償器）を用いて補うことアクティブ制御といいます。

演習A01…Flipped Classroom

$1^\circ$ 状態方程式 $\dot{x}(t)=ax(t)+bu(t)$ が与えられるときの漸近安定性の判定は、なぜ $\dot{x}(t)=ax(t)$ に基づいて行うのか説明せよ。

$2^\circ$ MATLABまたはSCILABに次のコマンドを与えて、どのグラフが漸近安定であるか述べよ。

MATLAB

%a01.m
clear all, close all
a1=1; a2=0; a3=-5; 
b=1; c=1;
sys1=ss(a1,b,c,[]);
sys2=ss(a2,b,c,[]);
sys3=ss(a3,b,c,[]);
t=0:0.01:1;
x0=1;
initial(sys1,sys2,sys3,x0,t)
grid
title('Stability of 1st-order System')
xlabel('time')
ylabel('x(t)')
legend('a=1','a=0','a=-5')
%eof

SCILAB

coming soon

図3 どのグラフが漸近安定であるか？

制御技術セミナー

投稿日時: 2021年10月7日投稿者: cacsd

現代制御　制御対象はダイナミックスの支配原理（物理）に逆らっては動けません。その支配方程式は微分方程式で記述されるので、制御対象のモデルもここを出発点とするのが自然です。制御とは「平衡状態の安定化」と捉えることができます。これらに沿う制御理論が現代制御（状態空間法）です。現代制御の理論にはストーリ性があって、漸近安定性、状態フィードバックと可制御性、状態オブザーバと可観測性、最適制御など、極めて美しく体系化されています。最適制御では、なんと未知変数が行列の場合の２次方程式（リッカチ方程式）を解きますが、現在の計算機を用いればこれは難なくこなせます。また、多入力・多出力系や閉ループ系の扱いもスマートです。これまで制御を意識したモノつくり（対象を敢えて不安定にして制御で安定化）が多くのイノベーションを生み出してきました。本セミナーでは現代制御の基礎知識を説明します。

ロバスト制御１（ゲインスケジューリング制御）　現代制御では制御対象のモデルを準備しなければなりませんが、一般には正確なモデルの構築は困難です。そこで、モデルの不確かさをどう扱うかが議論され、ポストモダンとしてのロバスト制御理論が発展してきました。特に、古典制御でのループ整形の考え方をベースにした混合感度問題への定式化がよく知られています。これはモデルの不確かさを入出力伝達特性の摂動によって間接的に扱う場合ですが、支配方程式のパラメータの不確かさを直接的に扱う場合もあります。パラメータが時間変動しない場合はラウスフルビッツの安定判別法が有用であり、変動パラメータの場合はゲインスケジューリングを行うことが試みられてきました。近年は線形行列不等式を制約条件とする最適化問題への定式化によって、混合感度問題やゲインスケジューリング問題の求解計算が容易になってきました。本セミナーではゲインスケジューリング制御の基礎知識を説明します。

ロバスト制御２（スライディングモード制御）　ロバスト制御に対するもう一つのアプローチがスライディングモード制御です。これは可変構造制御の一つで、線形制御則とスイッチング制御則からなります。瞬時瞬時に制御則を切り替えることにより、切り替えない場合に比べて優位性を生み出すことができます。スライディングモード制御は「スイッチング直線の近傍」まで状態を持っていく動作と、「スイッチング直線内に閉じ込める」２つの動作からなり、これによりモデルの不確かさに対応できるロバスト制御の特徴を持たせることができます。スイッチング関数の選択には現代制御のLQ設計法や固有ベクトル配置法が用いられ、リヤプノフ関数を構成することにより閉ループ系の安定性を保証することができます。本セミナーではスライディングモード制御の基礎知識を説明します。

感染制御

投稿日時: 2020年9月5日投稿者: cacsd

1. SIRモデルに基づく感染制御
1.1 SIRモデルの基礎式

感染ダイナミックスを表すSIRモデルでは、当該地域の人口 $N$ を

● $S(t)$ ：感受性(susceptible)状態の人口（未感染者数）
● $I(t)$ ：感染性(infectious)状態の人口（感染者数）
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口（回復者数）

に分け（ $N=S+I+R$ ）、これらが次の連立微分方程式に従うとします。

$\displaystyle{(1)\quad \left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right. }$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれます。

いま未感染者数 $S$ 、感染者数 $I$ 、回復者数 $R$ の当該地域の人口 $N$ に対する比率を考えると次式が成り立ちます。

$\displaystyle{(2)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1 }$

(1)と(2)から

$\displaystyle{(3)\quad \left\{\begin{array}{l} \underbrace{\frac{d\frac{I(t)}{N}}{dt}}_{\frac{di(t)}{dt}}=-\gamma \underbrace{\frac{I(t)}{N}}_{i(t)}+\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{I(t)}{N}}_{i(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{I}{N}}_{i(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{I(t)}{N}}_{i(t)} \end{array}\right. }$

すなわち、次の正規化（無次元化）されたSIRモデルを得ます。

$\displaystyle{(4)\quad \left\{\begin{array}{l} \frac{di(t)}{dt}=\gamma(-1+\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)}) i(t)\\ \frac{dr(t)}{dt}=\gamma i(t) \end{array}\right. }$

ここで、自粛率 $q(t)$ を導入し、 $\beta$ を $(1-q(t))\beta$ に置き換えて、次式を得ます。

$\displaystyle{(5)\quad \left\{\begin{array}{l} \frac{di(t)}{dt}=(-1+(1-q(t))\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)}) \gamma i(t)\\ \frac{dr(t)}{dt}=\gamma i(t) \end{array}\right. }$

図1にSIRモデル(5)に基づく感染シミュレーションの例を示します。ただし、自粛率 $q$ の時系列は８割削減と５割削減を適当に設定しています（ $1-q$ は削減率と呼ばれます）。

図1 正規化SIRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $q(t)=0,0.8,0.5$ ; $i(0)=0.003$ ）

1.2 SIRモデル（全国版）に基づく感染予測

図2に全国の新規感染者数の移動平均値を示します。これに基づく基本再生産数R0の値を図3に示します。これから7/25時点でのR0の値は $1.5$ 程度（ $\gamma=0.16$ ）と推測されます。

図2 全国の新規感染者数の実績値

図3 全国の基本再生産数の推定値

さらに感染予測結果を図4に示します。左図は7/1を初期時刻にとって、自粛は行われていないものとして、7/25現在の感染者数を予測したものです。右図は、7/25を初期時刻にとって、自粛率 $q$ を５割削減と２割削減を適当に設定して、向こう300日間の感染者数を予測したものです。新規感染者数の実績値は2000人を超えることはなかったので、もう少し厳しい自粛が行われていたと言えます。

図4 SIRモデルに基づく全国の感染予測の例

1.3 正規化SIRモデルに基づく感染制御

適当な累積感染者数と収束期間を実現するように自粛率 $q(t)$ の時系列を制御理論を用いて導出することを考えます。(5)の行列表示は次式となります。

$\displaystyle{(6)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t) \end{array}\right] =\gamma\left[\begin{array}{cc} 0 & 1\\ 0 & -1+(1-q(t))\bar{R}_0s(t) \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t) \end{array}\right] }$

さらに自粛率 $q(t)$ を１次遅れのダイナミックスを通して間接的に操作できるように政策変数（仮想的な操作変数） $u_q(t)$ を導入します。

$\displaystyle{(7)\quad \left\{\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t) \end{array}\right] =\left[\begin{array}{cc} 0 & \gamma\\ 0 & -\gamma+\gamma\bar{R}_0s(t) \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t) \end{array}\right] +\left[\begin{array}{c} 0\\ -\gamma\bar{R}_0s(t)i(t) \end{array}\right]q(t)\\ \dot{q}(t)=-\frac{1}{T_q}q(t)+\frac{1}{T_q}u_q(t) \end{array}\right. }$

すなわち、次式を得ます。

$\displaystyle{(8)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma+\gamma\bar{R}_0s(t)&-\gamma\bar{R}_0s(t)i(t)\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] +\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right]u_q(t) }$

これは２つの変動パラメータ $\alpha(t)=\bar{R}_0s(t)$ と $\beta(t)=\bar{R}_0s(t)i(t)$ をもつ、次のようなLPVモデル（線形パラメータ変動モデル）とみなすことができます。

$\displaystyle{(9)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] }_{\dot{x}(t)} =\left(\underbrace{\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma &0\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] }_{A_0} + \underbrace{\bar{R}_0s(t)}_{\alpha(t)} \underbrace{\left[\begin{array}{ccc} 0 & 0 & 0\\ 0 & \gamma & 0\\ 0 & 0 & 0 \end{array}\right] }_{A_1} \right.\\ \left.+ \underbrace{\bar{R}_0s(t)i(t)}_{\beta(t)} \underbrace{\left[\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & -\gamma\\ 0 & 0 & 0 \end{array}\right] }_{A_2} \right) \underbrace{\left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] }_{x(t)} +\underbrace{\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right] }_{B} \underbrace{u_q(t)}_{u(t)}\\ (0<\alpha(t),\beta(t)<\bar{R}_0) \end{array} }$

いま、２つの変動パラメータを適当な値に $\alpha^*=\bar{R}_0s^*$ と $\beta^*=\bar{R}_0s^*i^*$ に固定すると、次のLTIモデル（線形時不変モデル）を得ます。

$\displaystyle{(10)\quad \left[\begin{array}{c} \dot{r}(t)\\ \dot{i}(t)\\ \dot{q}(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma &0\\ 0 & -\gamma+\gamma\alpha^* &　-\gamma\beta^*\\ 0 & 0 &-\frac{1}{T_q} \end{array}\right] \left[\begin{array}{c} r(t)\\ i(t)\\ q(t) \end{array}\right] +\left[\begin{array}{c} 0\\ 0\\ \frac{1}{T_q} \end{array}\right]u_q(t) }$

そこで、２次形式評価関数

$\displaystyle{(11)\quad J=\int_0^\infty(r^2(t)+i^2(t)+q^2(t)+\rho^2u_q^2(t))dt }$

を最小化する状態フィードバック（LQ制御則）

$\displaystyle{(12)\quad u_q(f)=-f_1r(t)-f_2i(t)-f_3q(t) }$

の効果を確かめてみます。その結果を図5に示します。左図は自粛なしの場合、右図は時刻 $t=3$ 以降、LQ制御則を適用しています。ここで、 $s^*=0.85$ 、 $i^*=0.15$ 、 $T_q=0.2$ としています。時刻 $t=3$ ではすでに人口の15%（0.15 $N$ ）の感染者がいますので、強力な自粛が必要となり、 $q$ の値は１を超える大きな値となっています。そのためシミュレータ上では飽和要素を入れてこれを制限しています。結果的にいきなり10割削減相当、いわゆるロックダウンの自粛となっています。しかし最終的な累積感染者数は9割から３割程度に抑えられています。

図5 正規化SIRモデルに基づく感染制御（LQ制御）の例
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $i(0)=0.003$ ）

さて、 $\alpha_1\le\alpha(t)\le\alpha_2$ 、 $\beta_1\le\beta(t)\le\beta_2$ に対する次の内分式に注目します。

$\displaystyle{(13)\quad \begin{array}{l} \left[\begin{array}{c} \alpha(t)\\ \beta(t) \end{array}\right]= \underbrace{\frac{\alpha_2-\alpha(t)}{\alpha_2-\alpha_1}\frac{\beta_2-\beta(t)}{\beta_2-\beta_1}}_{p_{11}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_1\\ \beta_1 \end{array}\right]+ \underbrace{\frac{\alpha_2-\alpha(t)}{\alpha_2-\alpha_1}\frac{\beta(t)-\beta_1}{\beta_2-\beta_1}}_{p_{12}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_1\\ \beta_2 \end{array}\right]\\+ \underbrace{\frac{\alpha(t)-\alpha_1}{\alpha_2-\alpha_1}\frac{\beta_2-\beta(t)}{\beta_2-\beta_1}}_{p_{21}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_2\\ \beta_1 \end{array}\right]+ \underbrace{\frac{\alpha(t)-\alpha_1}{\alpha_2-\alpha_1}\frac{\beta(t)-\beta_1}{\beta_2-\beta_1}}_{p_{22}(\alpha(t),\beta(t))}\left[\begin{array}{c} \alpha_2\\ \beta_2 \end{array}\right] \\ \end{array} }$

$\displaystyle{(14)\quad p_{11}(\alpha(t),\beta(t))+p_{12}(\alpha(t),\beta(t))+p_{21}(\alpha(t),\beta(t))+p_{22}(\alpha(t),\beta(t))=1 }$

状態方程式(9)は、端点モデル

$\displaystyle{(15)\quad \begin{array}{l} \dot{x}(t)=\underbrace{A(\alpha_1,\beta_1)}_{A_{11}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_1,\beta_2)}_{A_{12}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_2,\beta_1)}_{A_{21}}x(t)+Bu(t)\\ \dot{x}(t)=\underbrace{A(\alpha_2,\beta_2)}_{A_{22}}x(t)+Bu(t) \end{array} }$

を、 $p_{11},p_{12},p_{12},p_{22}$ によって重み付けして、次式のように表されます。

$\displaystyle{(16)\quad \dot{x}(t)=\underbrace{(p_{11}(t)A_{11}+p_{12}(t)A_{12}+p_{21}(t)A_{21}+p_{22}(t)A_{22})}_{A(\alpha(t),\beta(t))}x(t)+Bu(t) }$

このとき４つの端点モデルを安定化する状態フィードバック

$\displaystyle{(17)\quad \begin{array}{l} u(t)=-F_{11}x(t)\\ u(t)=-F_{12}x(t)\\ u(t)=-F_{21}x(t)\\ u(t)=-F_{22}x(t) \end{array} }$

を同時安定化して求め（共通のリャプノフ関数を求め）、LPV制御則を次式で決めます。

$\displaystyle{(18)\quad u(t)=-\underbrace{(p_{11}(t)F_{11}+p_{12}(t)F_{12}+p_{21}(t)F_{21}+p_{22}(t)F_{22})}_{F(\alpha(t),\beta(t))}x(t) }$

この制御則の一例の効果を図6に示します。自粛率は１を超えることはありませんが、累積感染者数は９割から６割程度にしか抑制されていません。

図6 正規化SIRモデルに基づく感染制御（LPV制御）の例
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $i(0)=0.003$ ）

2. SIQRモデルに基づく感染制御
2.1 SIQRモデルの基礎式

SIRモデルは、次の３つの変数で表されていました。

● $S(t)$ ：感受性状態の人口（未感染者数）
● $I(t)$ ：感染性状態の人口（感染者数）
● $R(t)$ ：回復・隔離または免疫状態状態の人口（回復者数）

いま、感染者数を検査の有無によって、次のように区別します（ $I=H+Q$ ）。

● $H(t)$ ：未検査感染性人口（市中感染者数）
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

このとき、次式で表されるSIQRモデルを考えます。

$(101)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)H(t) \\ \frac{dH(t)}{dt}=-\gamma H(t)-pH(t)+(1-q)\beta S(t)H(t)\\ \frac{dQ(t)}{dt}=-\gamma Q(t)+pH(t)+q\beta S(t)H(t)\\ \frac{dR(t)}{dt}=\gamma (H(t)+Q(t)) \end{array}\right.$

ここで、SIQRモデルの第２式と第３式は、SIRモデルの第２式を、 $I=H+Q$ を用いて分けたものとなっています。その際、自粛率 $q$ の他に隔離率 $p$ も導入されていることに注意してください。 $pH(t)$ が日々発表される感染者数となります。

これら未感染者数、市中感染者数、隔離感染者数、回復者数の当該地域の人口 $N$ に対する比率を考えます。

$(102)\quad \underbrace{\frac{H(t)}{N}}_{h(t)}+\underbrace{\frac{Q(t)}{N}}_{q_u(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

このとき

$(103)\quad\left\{\begin{array}{l} \underbrace{\frac{d\frac{H(t)}{N}}{dt}}_{\frac{dh(t)}{dt}}=-\gamma \underbrace{\frac{H(t)}{N}}_{h(t)}-p \underbrace{\frac{H(t)}{N}}_{h(t)}+(1-q)\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{H(t)}{N}}_{h(t)}-\underbrace{\frac{Q(t)}{N}}_{q_u(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{H(t)}{N}}_{h(t)}\\ \underbrace{\frac{d\frac{Q(t)}{N}}{dt}}_{\frac{dq_u(t)}{dt}}=-\gamma \underbrace{\frac{Q(t)}{N}}_{q_u(t)}+p \underbrace{\frac{H(t)}{N}}_{h(t)}+q\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{H(t)}{N}}_{h(t)}-\underbrace{\frac{Q(t)}{N}}_{q_u(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{H(t)}{N}}_{h(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{H(t)}{N}}_{h(t)}+\gamma \underbrace{\frac{Q(t)}{N}}_{q_u(t)} \end{array}\right.$

すなわち、次の正規化（無次元化）されたSIQRモデルを得ます。

$(104)\quad\left\{\begin{array}{l} \frac{dh(t)}{dt}=-(\gamma+p)h(t)+(1-q)\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-h(t)-q_u(t)-r(t))}_{s(t)} \gamma h(t)\\ \frac{dq_u(t)}{dt}=-{\gamma}q_u(t)+{p}h(t)+q\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-h(t)-q_u(t)-r(t))}_{s(t)} \gamma h(t)\\ \frac{dr(t)}{dt}=\gamma(h(t)+q_u(t)) \end{array}\right.$

図7にSIQRモデル(104)に基づく感染シミュレーションの例を示します。ここで自粛率 $q$ の時系列は図1の半分にしています。これを隔離率 $p$ を導入することにより、図1と同程度の感染者数に抑え込むことができています。

図7 正規化SIQRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1; \bar{R}_0=2.5; q=0,0.8,0.5; p=0.1,0.5; i(0)=0.003$ ）

2.2 SIQRモデル（全国版）に基づく感染予測

図8 SIQRモデルに基づく全国の感染予測シミュレーション

2.3 正規化SIQRモデルに基づく感染制御

適当な累積感染者数と収束期間を実現するように自粛率 $q(t)$ と隔離率 $q(t)$ の時系列を制御理論を用いて導出してみます。(104)の行列表示は次式となります。

$(105)\quad\left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t) \end{array}\right] =\gamma\left[\begin{array}{ccc} 0 & 1 & 1\\ 0 & -1-\frac{p}{\gamma}+(1-q)\bar{R}_0s(t) & 0\\ 0 & \frac{p}{\gamma}+q\bar{R}_0s(t) & -1\\ \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t) \end{array}\right]$

自粛率 $q(t)$ と隔離率 $q(t)$ を１次遅れのダイナミックスを通して間接的に操作できるようにそれぞれの政策変数（仮想的な操作変数） $u_q(t)$ と $u_p(t)$ を導入します。

$(106)\quad\left\{\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t) \end{array}\right] =\left[\begin{array}{ccc} 0 & \gamma & \gamma\\ 0 & -\gamma+\gamma\bar{R}_0s(t) & 0\\ 0 & 0 & -\gamma\\ \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t) \end{array}\right]\\+ \left[\begin{array}{cc} 0 & 0\\ -\gamma\bar{R}_0s(t)h(t) & -1\\ \gamma\bar{R}_0s(t)h(t) & 1 \end{array}\right] \left[\begin{array}{c} q(t) \\ p(t) \end{array}\right]\\ \dot{q}(t)=-\frac{1}{T_q}q(t)+\frac{1}{T_q}u_q(t)\\ \dot{p}(t)=-\frac{1}{T_p}p(t)+\frac{1}{T_p}u_p(t) \end{array}\right.$

すなわち、次式を得ます。

$(107)\quad\begin{array}{l} \left[\begin{array}{c} \dot{r}(t)\\ \dot{h}(t)\\ \dot{q}_u(t)\\ \dot{q}(t)\\ \dot{p}(t) \end{array}\right] =\left[\begin{array}{ccccc} 0 & \gamma & \gamma&0 & 0\\ 0 & -\gamma+\gamma\bar{R}_0s(t) & 0&-\gamma\bar{R}_0s(t)h(t) & -1\\ 0 & 0 & -\gamma&\gamma\bar{R}_0s(t)h(t) & 1\\ 0&0&0&-\frac{1}{T_q}&0\\ 0&0&0&0&-\frac{1}{T_p} \end{array}\right] \left[\begin{array}{c} r(t)\\ h(t)\\ q_u(t)\\ q(t) \\ p(t) \end{array}\right]\\ +\left[\begin{array}{cc} 0&0\\ 0&0\\ 0&0\\ \frac{1}{T_q}&0\\ 0&\frac{1}{T_p} \end{array}\right] \left[\begin{array}{c} u_q(t) \\ u_p(t) \end{array}\right] \end{array}$

これは２つの変動パラメータ $\alpha(t)=\bar{R}_0s(t)$ と $\beta(t)=\bar{R}_0s(t)h(t)$ をもつ次のようなLPVモデル（線形パラメータ変動モデル）とみなすことができます。

PCR検査

投稿日時: 2020年9月5日投稿者: cacsd

[1] 現在、PCR検査拡大について賛否が分かれています。たとえば
●左翼メディアと煽り系専門家たちの罪
●無症状にPCR検査して隔離する世田谷モデルはどこがアホなのか説明する
●日本がコロナ第3波を避けられない決定的弱点
などを読むと、どれも説得力があって迷ってしまいます。ただ、一般的に言われている賛否の理由を整理してみると次のようになるかと思います。

●検査隔離拡大の反対理由：
(1) 偽陰性者が市中感染を広める（症状のある者・濃厚接触者に限定すべき）
(2) 偽陽性者数は偽陰性者数を上まわり、検査の無駄が多い（厚労省の特異度は99%、世界標準はほぼ100%）

●検査隔離拡大の賛成理由：
(3) 無症状者に繰り返し検査を行ない、検査陽性者を隔離することが、市中感染抑制につながる
(4) 無症状で検査陰性であることが経済を回すための国際ルールとなりつつある

個人的には(3)が重要だと思うのですが、今一つ理屈がよく分からないというのが、筆者ばかりでなく大方の率直な印象ではないでしょうか？その理解の一助となることを目指して、ここでは感染ダイナミックスに基づいて、一体どれくらいの人数を検査・隔離すればよいのか、具体的には市中感染者の何割ぐらいを検査・隔離すればよいのかを調べてみたいと思います。

[2] いま、当該地域の人口を $N$ 、感染率を $\beta$ 、回復率を $\gamma$ 、削減率を $q$ 、隔離率を $p$ とすると、感染が収束に向かうには、次式が成り立つ必要があります。

$\displaystyle{(5)\quad p>\gamma((1-q)R_0-1) }$

ここで、 $R_0=\frac{\beta N}{\gamma}$ は基本再生産数と呼ばれています。自粛をしていない場合は削減率 $q=0$ なので次式となります。

$\displaystyle{(6)\quad p>\gamma(R_0-1) }$

また市中感染者数を $H$ で表すと、これに隔離率 $p$ をかけた $p\times H$ が日々報告される感染者数となります。したがって、日々報告される感染者数を $p$ で割れば、検査すべき市中感染者数 $H$ を見積もることができます。

いま(5)の右辺を

$\displaystyle{(7)\quad p^*(R_0,q,\gamma)=\gamma((1-q)R_0-1) }$

とおいて、 $R_0$ を1.5, 2.0, 2.5, 3.0の４通り、 $q$ を0, 0.2, 0.5, 0.8の４通り、 $\gamma$ を0.04, 0.08, 0.12, 0.16の４通り、すべての組合せ $4^3=64$ 通りについて、 $p^*(R_0,q,\gamma)$ の値を計算したものを、表1の左側に示します（負値を与える組合せはすでに収束に向かっていることを意味します）。また対応する $1/p^*(R_0,q,\gamma)$ を表1の左側に示します。これは検査拡大率を意味します。その棒グラフを描いた図1から次が分かります。

●基本再生数 $R_0$ が小さいほど検査拡大率は大きい
●削減率 $q$ が大きいほど検査拡大率は大きい

したがって、検査拡大は、基本再生数 $R_0$ が小さい（感染力が小さい）ときほど、削減率 $q$ が大きい（自粛レベルが高い）ときほど必要と言えます。すなわちよく言われるように感染が一旦収まったときに徹底した検査を行うべきという主張を裏付けるものと思われます。

表1 $p^*(R_0,q,\gamma)$ と $1/p^*(R_0,q,\gamma)$

図1 表1の $1/p^*(R_0,q,\gamma)$ の棒グラフ

注感染ダイナミックスの著名なものはSIRモデルと言われます。これに自粛率 $1-q$ （ $q$ は削減率）ばかりでなく、隔離率 $p$ を導入し、SIQRモデルとして提案されたのは小田垣先生です。詳しくは感染ダイナミックスをご覧ください。そこで示している初期成長率(31)式から上の不等式(5)が得られます。あくまで感染初期に適用できる式です。

[3] さて、図1に全国の新規感染者数の実績値の時系列を示します。これを見ると、５月下旬から６月までは第１波が収束していた時期です。ところが７月から急速に感染者数が増加していきます。７月は正に自粛なしの期間と言ってよいかと思います。図2に対応する基本生産数を図3に示します。７月の１か月間では、 $R_0=1.5$ 、 $\gamma=0.12$ と読み取ることができ、(7)は次式となります。

$\displaystyle{(7)\quad p^*(1.5,0,0.12)=\gamma(R_0-1)=0.12\times 0.5=0.06 }$

したがって、全国レベルで日々1000人前後の感染者数が報告されていたので、この約17倍(1/0.06=17)の検査が毎日必要だったと言えそうです。さらに３か月間週１回のペースで繰り返し12回検査するとすると、約20万人規模の検査が必要だったと言えそうです。

以上は、少し大雑把で強引な議論ですが、経済を回すために自粛せずに、感染収束に向かわせるためにはかなりの検査規模が必要になることが分かります。このことが、反対理由(1),(2)から生じるデメリット（検査の無駄）を打ち消すことになるのかは明確ではありませんが、ご参考になればと思います。

図3 全国の新規感染者数の実績値

図4 全国の基本再生産数の推定値

感染ダイナミックス（2020/7/28)

投稿日時: 2020年7月28日投稿者: cacsd

感染ダイナミックスはSIRモデルによって表すことができます。図1は、これをMATLAB/Simulinkで表したものです。

図1 正規化SIRモデルに基づく感染シミュレータの例(MATLAB/Simulink)
（ $\gamma=1$ ; $\bar{R}_0=2.5$ ; $q(t)=0,0.8,0.5$ ; $i(0)=0.003$ ）

図2に全国の新規感染者数の移動平均値を示します。これに基づく基本再生産数 $\bar{R}_0$ の値を図3に示します。これから7/25時点での $\bar{R}_0$ の値は $1.5$ 程度（ $\gamma=0.16$ ）と推測されます。さらに予測結果を図4に示します。左図は7/1を初期時刻にとって、自粛は行われていないものとして、7/25現在の感染者数を予測したものです。右図は、7/25を初期時刻にとって、自粛率 $q$ を５割削減と２割削減を適当に設定して、向こう300日間の感染者数を予測したものです。

図2 全国の新規感染者数の実績値

図3 全国の基本再生産数の推定値

図4 SIRモデルに基づく全国の感染予測の例

これから早急に自粛要請をかけないと感染者数は爆発的に増加していくことが分かります。ただ自粛率という数値を行動変容にどう結び付けるかは不明です。我々の感染拡大への想像力と感染対策の知恵が問われています。

感染ダイナミックス（2020/6/26）

投稿日時: 2020年6月26日投稿者: cacsd

ここでは次の文献にもとづいて「感染ダイナミックス」について述べます。
[1] 稲葉寿編著：「感染症の数理モデル」,培風館（2008）
[2] 牧野淳一郎：「3.11以後の科学リテラシーno.89」,科学, vol.90, no.5, 岩波書店(2020)
[3] 小田垣孝：「新型コロナウイルスの蔓延に関する一考察」
なお、感染ダイナミックスの数理モデルは、次の文献が先駆けとなっています。
Kermack, W. O. and McKendrick, A. G.（1927）. Contributions to the mathematical theory of epidemics – I, Proceedings of the Royal Society Series A, 115, 700–721

1 SIRモデル
1.1 基本再生産数の概念
感染ダイナミックスの数理モデルのうち最も基本的なものはSIRモデルとして知られています。ここでは文献[1]に従って、その概要を述べます。

いま当該地域の人口 $N$ は次のような構成であるとします。

$(1)\quad S(t)+I(t)+R(t)=N$

ただし、時刻 $t$ において

● $S(t)$ ：感受性(susceptible)状態の人口
● $I(t)$ ：感染性(infectious)状態の人口
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口

以下では $S$ 、 $I$ 、 $R$ をそれぞれ未感染者数、感染者数、回復者数と呼びます。いくつかの仮定のもと、感染ダイナミックスの支配方程式として、次のSIRモデルを考えます。

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれています。

いま時刻 $t$ における感染状況を３項組 $(S(t),I(t),R(t))$ で表します。感染前の自明な平衡状態は

$(3)\quad (S,I,R)=(N,0,0)$

ですが、これが初期時刻 $t=0$ において、感染者がごく少数 $I_0$ だけ現れ

$(4)\quad (S_0,I_0,R_0)=(N-I_0,I_0,0)$

となったとします。このときの関心事は、感染流行が始まるか、それとも感染収束するかどうかです。これは(2)の第２式において、感染初期では $S(t)\simeq S_0$ として

$(5)\quad \frac{dI(t)}{dt}=\underbrace{(\beta S_0-\gamma)}_{\lambda_0}I(t), \quad I(0)=I_0$

の解

$(6)\quad I(t)=\exp(\lambda_0 t)I_0$

の振舞いを調べることになります。明らかに、初期成長率 $\lambda_0$ の正負によって、感染流行か感染収束かが決まります。これを

$(7)\quad \lambda_0=\beta S_0-\gamma=\gamma(\frac{\beta S_0}{\gamma}-1)$

と変形して、ここで現れた定数

$(8)\quad R_0=\frac{\beta S_0}{\gamma}$

を基本再生産数と名付け

● $R_0>1$ のとき感染流行
● $R_0<1$ のとき感染収束

となると判断します。この方法は閾値原理と呼ばれています。

さて、基本再生産数は次のように表されます。

$(9)\quad R_0=\beta S_0\underbrace{\frac{1}{\gamma}}_{T}=1+\underbrace{(\beta S_0-\gamma)}_{\lambda_0}\underbrace{\frac{1}{\gamma}}_{T}$

ここで、 $T$ は感染性時間（感染力をもつ時間）を表すと解釈できます。実際、感染力がなくなったと仮定すると( $\beta=0$ )、(2)の第２式は

$(10)\quad\frac{dI(t)}{dt}=-\gamma I(t), \quad I(0)=I_0$

となり、この解は

$(11)\quad I(t)=\exp(-\gamma t)I_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、 $T=1/\gamma$ の5～10倍の時間が経てば、感染者数 $I_0$ は十分零になるからです。

そこで(9)において、 $\beta S_0$ は１人が単位時間に２次感染させる人数、これに感染性時間 $T$ をかけるので、基本再生産数は感染者１人が２次感染させる人数を表すと言えます。

その算出には初期成長率 $\lambda_0$ が必要ですが、次式を用います。

$(12)\quad \frac{I(t+\Delta t)}{I(t)}=\frac{\exp(\lambda_0 (t+\Delta t))I_0}{\exp(\lambda_0 t)I_0}=\exp(\lambda_0 \Delta t) \Rightarrow \lambda_0=\ln(\frac{I(t+\Delta t)}{I(t)})/\Delta t$

1.2 基本再生産数の推定例
専門家による基本再生産数の推定法は、感染学上の知見に基づく統計推測手法によるものと思われますが、詳細は明らかではありません。そこで、次のような基本再生産数の推定手法を考えてみます。

基本再生産数 $R_0$ の推定手法１：
$1^\circ$ 当該地域の感染者数データの移動平均をとる
$2^\circ$ 次式を用いて時刻 $t$ における（初期）成長率 $\lambda_0$ を計算する

$(13)\quad\lambda_0(t)=\ln(\frac{I(t)}{I(t-\Delta t)})/\Delta t$

$3^\circ$ 次式を用いて時刻 $t$ における基本再生産数 $R_0$ を計算する

$(14)\quad R_0(t)=1+\lambda_0(t) T$

図1に全国、東京都、福岡県、NYの感染数を7日間移動平均したデータを示します。図2に日ごとに基本再生産数 $R_0$ を推定した例を示します。ただし、 $2^\circ$ で $\Delta t=7$ としています。 $3^\circ$ で回復率は $\gamma=0.04,0.08,0.12,0.16$ の４通りを試し、 $R_0<0$ となる場合は0としています。

図1 感染者数の実データ（移動平均後）
国内データソース：NHK 特設サイト新型コロナウイルス
NYデータソ－ス：New York State Department of Health COVID-19 Tracker

図2 基本再生産数 $R_0$ の推定例

図2の妥当性を検討するために、「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日）に記載されている次の公表値に注目します。

●全国の実効再生産数：　 $R_0(3/25)=2$ , 　 $R_0(4/10)=0.7$
●東京都の実効再生産数： $R_0(3/14)=2.6$ , $R_0(4/10)=0.5$

ここで基本再生産数の代わりに実効再⽣産数という術語が用いられています。ここでは両者の理論上の違いまでには踏み込まず、実効再⽣産数は基本再生産数 $R_0$ に実質的に等しいとして話を進めます。上記の値は上記資料に掲載されている次の図版から読み取ることができます。

図3 専門家会議資料における実効再生産数の推定

図3の左側は全国について、右側は東京都について、感染者数と実効再生産数が示されています。下側に示されている棒グラフは上側の感染者数の移動平均をとったものと思われ、これに基づいて実効再生産数が計算されています。ただし、移動平均後のデータは２週間ほど前倒しされています（ピーク値と $R_0$ が１を切る時点が少しずれています）。実際、全国（東京都）の感染者数のピークは4/11（4/9）ですが、移動平均後は3/28（3/27）にピークが現れています。図2において全国の基本再生産数が0.7まで下がるのは早くて4/20当たり、東京都の基本再生産数が0.5まで下がるのは早くて4/23当たりですから、ほぼ２週間遅れとなっています。また、回復率 $\gamma$ の値に依らず、感染者数がピークを取る時点で $R_0$ が１を切り、収束に向かっている傾向が捉えられています。

問題はどの回復率 $\gamma$ が妥当かですが、文献でよく見るのは $\gamma=0.04$ です。しかしこれでは直近で $R_0$ が負となりました（ $\gamma=0.12$ , $0.16$ では常に正でした）。回復率の逆数は感染性時間 $T=1/\gamma$ ですから、 $\gamma=0.04$ , $0.08$ , $0.12$ , $0.16$ に対応して、 $T=25$ , $12.5$ , $8.3$ , $6.25$ です。指数関数的に減衰する場合は、初期値の10%まで減衰する期間は $2.3T$ 、2%まで減衰する期間は $3.9T$ ですから、たとえば、 $\gamma=0.16$ , $T=6.25$ でもそうおかしくはないかと思います。

2 感染シミュレーション
2.1 無次元化SIRモデル
前章では、SIRモデルのうち第２式のみを用いて議論しました。ここでは３本の微分法方程式を連立させて、 $S$ 、 $I$ 、 $R$ の時系列のシミュレーションを行います。

$S$ 、 $I$ 、 $R$ の間には(1)の拘束があるので、実質的には２本の微分方程式を連立させて解きます。(2)の第２式と(1)からSを消去して、第３式を合わせると次式を得ます（第１式はこれらから導出できます）。

$(15)\quad\left\{\begin{array}{l} \frac{dI(t)}{dt}=-\gamma I(t)+\beta \underbrace{(N-I(t)-R(t))}_{S(t)}I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで３つのパラメータ $N$ 、 $\beta$ 、 $\gamma$ がありますが、無次元化という操作を行って、１つのパラメータにまとめます。これは文献[2]を参考にしています。

まず、人口 $N$ への依存を除くために、各変数の割合（人口比）、すなわち、感染者数の割合、回復者数の割合、未感染者数の割合を、それぞれ

$(16)\quad i(t)=I(t)/N, r(t)=R(t)/N, s(t)=S(t)/N$

のように小文字で表します。(1)より次式が成り立ちます。

$(17)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

これを用いて、(15)から次式を得ます。

$(18)\quad\left\{\begin{array}{l} \underbrace{\frac{d\frac{I(t)}{N}}{dt}}_{\frac{di(t)}{dt}}=-\gamma \underbrace{\frac{I(t)}{N}}_{i(t)}+\beta N(\underbrace{\frac{N}{N}}_{1}-\underbrace{\frac{I(t)}{N}}_{i(t)}-\underbrace{\frac{R(t)}{N}}_{r(t)})\underbrace{\frac{I}{N}}_{i(t)}\\ \underbrace{\frac{d\frac{R(t)}{N}}{dt}}_{\frac{dr(t)}{dt}}=\gamma \underbrace{\frac{I(t)}{N}}_{i(t)} \end{array}\right.$

次に、感染性時間 $T=1/\gamma$ が単位時間となるように、次のように変形します。

$(19)\quad\left\{\begin{array}{l} \underbrace{\frac{di(t)}{d\frac{t}{T}}}_{\frac{di(t)}{d\tau}}=-i(t)+\underbrace{\frac{\beta N}{\gamma}}_{\bar{R}_0} \underbrace{(1-i(t)-r(t))}_{s(t)} i(t)\\ \underbrace{\frac{dr(t)}{d\frac{t}{T}}}_{\frac{dr(t)}{d\tau}}=i(t) \end{array}\right.$

すなわち、 $\tau=t/T$ とおいて、次式を得ます。

$(20)\quad\left\{\begin{array}{l} {\frac{di(\tau)}{d\tau}}=-i(\tau)+\bar{R}_0 s(\tau)i(\tau)=\underbrace{(-1+\bar{R}_0 s(\tau))}_{\lambda(\tau)}i(\tau)\\ {\frac{dr(\tau)}{d\tau}}=i(\tau) \end{array}\right.$

ここで、次の定数が定義されていることに注意してください。

$(21)\quad \bar{R}_0=\frac{\beta N}{\gamma}>\frac{\beta S_0}{\gamma}$

すなわち、 $\bar{R}_0$ は基本再生産数 $R_0$ の上界を与えています。ただし、感染初期の実データでは $S_0\simeq N$ が成り立ち、 $\bar{R}_0\simeq R_0$ です。以下では、 $R_0$ は $\bar{R}_0$ の意味で使用しています。

注1：(20)における $i(\tau)$ , $s(\tau)$ , $r(\tau)$ , $\lambda(t)$ はそれぞれ $i(\tau T)$ , $s(\tau T)$ , $r(\tau T)$ , $\lambda(\tau T)$ として、離散時間の差分方程式として表すべきかもしれませんが、以下では上式を用います。

2-2 感染シミュレーション
SIRモデル(20)の解の振舞いは、数値計算を行って容易にチェックできます。図3は $i(0)=0.001$ のとき、基本再生産数 $R_0$ の値によって感染ダイナミックスがどのように変わるかを示しています。青線が未感染者数 $S$ 、赤線が感染者数 $I$ 、緑線が回復者数（＝累積感染者数） $R$ の割合です。これから、基本再生産数 $R_0$ の値によって、最終的にどれくらいの割合で感染者が出るか、どれくらいの収束期間がかかるかを読み取ることができます。ただし、どの図も $R_0$ が一定の場合であることに注意してください。これらの図から、もし $R_0$ の値を順次下げていくことができれば、最終的な感染者数の割合は小さくなることが分かります。

図4 基本再生産数 $R_0$ に依存する感染ダイナミックス

注2：図4において、時間軸の単位は感染性時間 $T=1/\gamma$ です。したがって、収束期間を見積もるためには、回復率 $\gamma$ の推定が必要になります。

図5は、いま基本再生産数は $R_0=2.5$ であるとし、ある時点からこれを $(1-q)R_0$ に低減したとき、感染者数がどのように変化するかを示しています( $i(0)=0.001$ の場合です)。 $q=0.8$ のときが8割低減、 $q=0.2$ のときが2割低減です。

図5 基本再生産数 $R_0=2.5$ を $(1-q)R_0$ に低減したときの効果

注3：「８割削減の妥当性」についてですが、これはツイート_R020405が基礎となっているようです。ここでは基本再生産数 $R_0=2.5$ の場合であることが表明されていますので図5が参考になります。これは「基本再生産数を８割減することが必要で、２割減では効果がない」という主張をほぼ裏付けていると言えます。

2-2 感染初期における予測
実データによる感染予測ではSIRモデル(20)を解くわけですが、その際 $R_0$ を推定する必要があります。そのために次の問題設定を考えます。

基本再生産数 $R_0$ の推定手法２：
いま人口を $N$ 、感染者数の時系列を $I(t_1),\cdots,I(t_n)$ とします。またパラメータ $R_0$ を含むSIRモデル(20)の解の $t_1,\cdots,t_n$ における値を $i(t_1),\cdots,i(t_n)$ とします。このとき評価関数

$(22)\quad J(R_0)=(I(t_1)/N-i(t_1))^2+\cdots+(I(t_n)/N-i(t_n))^2$

を最小化するようなパラメータ $R_0$ の値 $R_0^*$ を、数値最適化の手法を用いて求めます。

この手法を用いて、全国の1/24～5/5の新規感染者数の一部のデータ（1/24～4/7、赤色〇）をもとに、基本再生産数を推定したところ $R_0^*$ =1.08となりました。図6は、感染者数の実績値（〇）と、SIRモデルによる計算値（青色実線）を重ねてプロットしたものです。感染初期の感染者数の指数関数的増加傾向がよく捉えられています。
図6 日本における感染者数の推移

首相は4/7日夜の記者会見で『このペースで感染拡大が続けば（感染者が）２週間後には１万人、１カ月後には８万人を越える』と指摘しました。図7は図6の感染者数の累積値を示しています。起点となる日を4/7（75日目）としますと、２週間後（88日目）より早く１万人を超えています（実際に１万人を超えたのは4/18（86日目）でした）。また８万人に到達するのは112日目あたりで、１か月後(105日目)より少し遅くなっています。
図7 日本における累積感染者数の推移

注4：上の問題を解くとき、少数の感染者が出た日を第１日として、 $t_1=1$ 、以後 $t_2=2,t_3=3\cdots$ とします。一方、SIRモデル(20)の解も $t_1=1,t_2=2,\cdots$ において求めます。この方法で感染初期の指数関数的増大の様子をよく捉えることができます。一方、4/11放映のNHKスペシャルより、基本再生産数 $R_0$ は大体1.5程度であることが表明されています。上で求めた $R_0^*$ の値1.08は公表値1.5より低めに出ています。そこで、日ごとの実データに回帰をとる場合 $T=1$ と考えられるので、(9)より $初期成長率\lambda_0$ を求めると

$(23)\quad \underbrace{1.08}_{R_0^*}=1+\lambda_0\times\underbrace{1}_{T}\quad\Rightarrow\quad \lambda_0=0.08}$

となります。これをたとえば感染性時間 $T=1/\gamma=1/0.16=6.25$ を用いて

$(24)\quad R_0=1+\underbrace{0.08}_{\lambda_0}\times\underbrace{6.25}_{T}=1.5}$

のように修正すると、公表値と一致します。

3. SIQRモデル
5/5に小田垣先生が新型コロナウイルスの蔓延に関する一考察を発表されました。現在指摘されている「PCR検査不足」の主張の妥当性を議論する際のベースを与えることのできるきわめて興味深い論文です。以下では、本稿の枠組みの中でフォローアップさせていただきたいと思います。

当該論文では、感染性人口 $I(t)$ を

● $H(t)$ ：未検査感染性人口（市中感染者数）　⇒　筆者の方で記号をIからHに変更
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

に分けて、次のSIQRモデルが提案されています。

$(25)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)H(t) \\ \frac{dH(t)}{dt}=-\gamma H(t)-pH(t)+(1-q)\beta S(t)H(t)\\ \frac{dQ(t)}{dt}=-\gamma Q(t)+pH(t)+q\beta S(t)H(t)\\ \frac{dR(t)}{dt}=\gamma (H(t)+Q(t)) \end{array}\right.$

これをSIRモデル

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

と比較します。それぞれの第 $i$ 式を、SIRi、SIQRiと参照します。

●SIQR4はSIR3に $I=H+Q$ を代入しただけです。
●SIQR1はSIR1の $I=H+Q$ のうち隔離感染者数 $Q$ の影響はないとしたものです。
●SIQR2とSIQR3は、SIR2の $I=H+Q$ を分けたものに相当します。SIR2の $\beta S(t)I(t)$ については

$(26)\quad \beta S(t)I(t)=(1-q)\beta S(t)I(t)+q\beta S(t)I(t)$

と分解し、 $I=H+Q$ のうち隔離感染者数 $Q$ の影響はないとしています。ここでパラメータ $q$ は感染率 $\beta$ に対する削減率とみなせます。またSIQR2とSIQR3には、 $pH(t)$ が差し引きしてあります。この $pH(t)$ は検査後陽性判定を受け隔離される人口で、日々発表される感染者数とされています。ここでパラメータ $p$ は市中感染者に対する隔離率とみなせます。

このとき、SIRモデルの場合の(5)と同様に

$(27)\quad \frac{dH(t)}{dt}=\underbrace{((1-q)\beta S_0-\gamma-p)}_{\lambda_0(q,p)}H(t), \quad H(0)=H_0$

が成り立ち、初期成長率 $\lambda_0(q,p)$ は

$(28)\quad \lambda_0(q,p)=(\gamma+p)((1-q)\underbrace{\frac{\beta S_0}{\gamma+p}}_{\tilde{R}_0}-1)$

と表すことができます。これをSIRモデルの場合の

$(7)\quad \lambda_0=\gamma(\underbrace{\frac{\beta S_0}{\gamma}}_{R_0}-1)$

と比較してみますと、 $\gamma$ が $\gamma+p$ に代わっており、基本再生産数が小さくなっています。したがって、SIQRモデルによれば、隔離率 $p$ を上げれば基本再生産数が小さくなる可能性を理解できます。

(27)の解は

$(29)\quad H(t)=\exp(\lambda_0(q,p)t)H_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、初期値 $H_0$ の10%まで下がる時間は( $\lambda_0(q,p)<0$ に注意して)

$(30)\quad \exp(\lambda_0(q,p)t)=0.1\quad\Rightarrow\quad t=\frac{\ln(0.1)}{|\lambda_0(q,p)|}=\frac{2.3}{|\lambda_0(q,p)|}$

で求められます。

当該論文の重要な貢献は、初期成長率を

$(31)\quad \lambda_0=\beta S_0-\gamma\quad\Rightarrow\quad \lambda_0(q,p)=(1-q)\beta S_0-(\gamma+p)$

のように、隔離率 $p$ と削減率 $q$ で変化させることができることを明らかにしていることです。そして、当該論文では、全国の感染者数データについて、次の結論を示しています。

Case 1) 現状： $\lambda_0(0.6,p)=-0.075$ 　⇒　収束期間31日
Case 2) ８割自粛： $\lambda_0(0.8,p)=-0.101$ 　⇒　収束期間23日
Case 3) ロックダウン： $\lambda_0(1,p)=-0.126$ 　⇒　収束期間18日
Case 4) 検査４倍増： $\lambda_0(0,4p)=-0.288$ 　⇒　収束期間8日
Case 5) ５割自粛＋検査２倍増： $\lambda_0(0.5,2p)=-0.159$ 　⇒　収束期間14日

これらは次を(31)と(30)に代入して確認できます。

$(32)\quad \beta S_0=0.126$ 、 $\gamma=0.03$ 、 $p=0.096$

たとえば、Case 1)については

$(33)\quad \underbrace{0.4}_{1-q}\times \underbrace{0.126}_{\beta S_0}-\underbrace{0.096}_{p}-\underbrace{0.03}_{\gamma}=\underbrace{-0.075}_{\lambda_0}\quad\Rightarrow\quad \underbrace{2.3/0.075}_{2.3/|\lambda_0|}=31$

当該論文では、基本再生産数を

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.126}{0.03+0.096}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.051}{0.03+0.096}=0.405$

のように計算しています。これから現状の削減率は $q=1-0.051/0.126=0.595$ であることが確かめられます。また、回復率としては $\gamma=0.03$ 、感染性期間が $T=1/\gamma=33$ が採用されています。さらに、感染減少期間に入る前の感染頂上期間では $\tilde{R}_0=1$ となるように、隔離率 $p=0.096$ 、 $1/p=10.4$ が定められています。このように、感染者数データから隔離率 $p$ を見積もることができ、 $1/p$ 倍して市中感染者数を推定できるという意味で大変興味深いと思います。

注5：本稿の計算では、基本再生産数は

感染増大期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.128}{0.04}=3.191$
感染減少期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.022}{0.04}=0.545$

となります。そこで当該論文と同様の計算をしてみますと

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.128}{0.04+0.088}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.022}{0.04+0.088}=0.545$

となります。これから削減率 $q=1-0.022/0.128=0.829$ 、隔離率 $p=0.088$ 、 $1/p=11.4$ が分かります。これらを用いると次の計算結果を得ます。

Case 2) ６割自粛： $\lambda_0(0.6,p)=-0.077$ 　⇒　収束期間30日
Case 1) 現状： $\lambda_0(0.829,p)=-0.106$ 　⇒　収束期間22日
Case 3) ロックダウン： $\lambda_0(1,p)=-0.128$ 　⇒　収束期間18日
Case 4) 検査４倍増： $\lambda_0(0,4p)=-0.263$ 　⇒　収束期間9日
Case 5) ５割自粛＋検査２倍増： $\lambda_0(0.5,2p)=-0.151$ 　⇒　収束期間15日

本稿のデータ分析では、すでに８割強の自粛が行われている結果になっていますが、PCR検査を増やすことの効果を確認できたと言えます。ただし、偽陽性と偽陰性の可能性を勘案して、陽性の出る可能性の高い場合に検査を行う必要があることは述べるまでもありません。

感染ダイナミックス（2020/5/17)

投稿日時: 2020年5月17日投稿者: cacsd

このサイトはコロナ禍が始まって最初にまとめました。ただし、図1と図2は随時更新しています。また感染制御についてもまとめようとしています。

1 はじめに
新型コロナウィルスに関する検討結果を公開しています。「なぜ外出自粛が必要かを説得力をもって説明することは難しい」という問題意識から検討を始めました（ツイート_R020330、ツイート_R020402）。また「８割削減」に示される数値の根拠に感心が集まっています。筆者の専門は感染学ではなく、システム制御工学ですが、一定の説明はできると思っています。

ポイントは「感染ダイナミックス」への理解です。ここでは次の文献を参照します。
[1] 稲葉寿編著：「感染症の数理モデル」,培風館（2008）
[2] 牧野淳一郎：「3.11以後の科学リテラシーno.89」,科学, vol.90, no.5, 岩波書店(2020)
なお、感染ダイナミックスの数理モデルは、次の文献が先駆けとなっています。
Kermack, W. O. and McKendrick, A. G.（1927）. Contributions to the mathematical theory of epidemics – I, Proceedings of the Royal Society Series A, 115, 700–721

2 SIRモデル
2.1 基本再生産数の概念
感染ダイナミックスの数理モデルのうち最も基本的なものはSIRモデルとして知られています。ここでは文献[1]に従って、その概要を述べます。

いま当該地域の人口 $N$ は次のような構成であるとします。

$(1)\quad S(t)+I(t)+R(t)=N$

ただし、時刻 $t$ において

● $S(t)$ ：感受性(susceptible)状態の人口
● $I(t)$ ：感染性(infectious)状態の人口
● $R(t)$ ：回復・隔離または免疫状態(recovered/removed/immune)状態の人口

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

ここで、 $\beta$ は感染率、 $\gamma$ は回復率と呼ばれています。

いま時刻 $t$ における感染状況を３項組 $(S(t),I(t),R(t))$ で表します。感染前の自明な平衡状態は

$(3)\quad (S,I,R)=(N,0,0)$

ですが、これが初期時刻 $t=0$ において、感染者がごく少数 $I_0$ だけ現れ

$(4)\quad (S_0,I_0,R_0)=(N-I_0,I_0,0)$

$(5)\quad \frac{dI(t)}{dt}=\underbrace{(\beta S_0-\gamma)}_{\lambda_0}I(t), \quad I(0)=I_0$

の解

$(6)\quad I(t)=\exp(\lambda_0 t)I_0$

の振舞いを調べることになります。明らかに、初期成長率 $\lambda_0$ の正負によって、感染流行か感染収束かが決まります。これを

$(7)\quad \lambda_0=\beta S_0-\gamma=\gamma(\frac{\beta S_0}{\gamma}-1)$

と変形して、ここで現れた定数

$(8)\quad R_0=\frac{\beta S_0}{\gamma}$

を基本再生産数と名付け

● $R_0>1$ のとき感染流行
● $R_0<1$ のとき感染収束

となると判断します。この方法は閾値原理と呼ばれています。

さて、基本再生産数は次のように表されます。

$(9)\quad R_0=\beta S_0\underbrace{\frac{1}{\gamma}}_{T}=1+\underbrace{(\beta S_0-\gamma)}_{\lambda_0}\underbrace{\frac{1}{\gamma}}_{T}$

ここで、 $T$ は感染性時間（感染力をもつ時間）を表すと解釈できます。実際、感染力がなくなったと仮定すると( $\beta=0$ )、(2)の第２式は

$(10)\quad\frac{dI(t)}{dt}=-\gamma I(t), \quad I(0)=I_0$

となり、この解は

$(11)\quad I(t)=\exp(-\gamma t)I_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、 $T=1/\gamma$ の5～10倍の時間が経てば、感染者数 $I_0$ は十分零になるからです。

その算出には初期成長率 $\lambda_0$ が必要ですが、次式を用います。

$(12)\quad \frac{I(t+\Delta t)}{I(t)}=\frac{\exp(\lambda_0 (t+\Delta t))I_0}{\exp(\lambda_0 t)I_0}=\exp(\lambda_0 \Delta t) \Rightarrow \lambda_0=\ln(\frac{I(t+\Delta t)}{I(t)})/\Delta t$

2.2 基本再生産数の推定例
専門家による基本再生産数の推定法は、感染学上の知見に基づく統計推測手法によるものと思われますが、詳細は明らかではありません。そこで、次のような基本再生産数の推定手法を考えてみます。

$(13)\quad\lambda_0(t)=\ln(\frac{I(t)}{I(t-\Delta t)})/\Delta t$

$3^\circ$ 次式を用いて時刻 $t$ における基本再生産数 $R_0$ を計算する

$(14)\quad R_0(t)=1+\lambda_0(t) T$

図2 基本再生産数 $R_0$ の推定例

図2の妥当性を検討するために、「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日）に記載されている次の公表値に注目します。

●全国の実効再生産数：　 $R_0(3/25)=2$ , 　 $R_0(4/10)=0.7$
●東京都の実効再生産数： $R_0(3/14)=2.6$ , $R_0(4/10)=0.5$

図3 専門家会議資料における実効再生産数の推定

3 感染シミュレーション
3.1 無次元化SIRモデル
前章では、SIRモデルのうち第２式のみを用いて議論しました。ここでは３本の微分法方程式を連立させて、 $S$ 、 $I$ 、 $R$ の時系列のシミュレーションを行います。

$(15)\quad\left\{\begin{array}{l} \frac{dI(t)}{dt}=-\gamma I(t)+\beta \underbrace{(N-I(t)-R(t))}_{S(t)}I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

まず、人口 $N$ への依存を除くために、各変数の割合（人口比）、すなわち、感染者数の割合、回復者数の割合、未感染者数の割合を、それぞれ

$(16)\quad i(t)=I(t)/N, r(t)=R(t)/N, s(t)=S(t)/N$

のように小文字で表します。(1)より次式が成り立ちます。

$(17)\quad \underbrace{\frac{I(t)}{N}}_{i(t)}+\underbrace{\frac{R(t)}{N}}_{r(t)}+\underbrace{\frac{S(t)}{N}}_{s(t)}=1$

これを用いて、(15)から次式を得ます。

次に、感染性時間 $T=1/\gamma$ が単位時間となるように、次のように変形します。

すなわち、 $\tau=t/T$ とおいて、次式を得ます。

ここで、次の定数が定義されていることに注意してください。

$(21)\quad \bar{R}_0=\frac{\beta N}{\gamma}>\frac{\beta S_0}{\gamma}$

3-2 感染シミュレーション
SIRモデル(20)の解の振舞いは、数値計算を行って容易にチェックできます。図3は $i(0)=0.001$ のとき、基本再生産数 $R_0$ の値によって感染ダイナミックスがどのように変わるかを示しています。青線が未感染者数 $S$ 、赤線が感染者数 $I$ 、緑線が回復者数（＝累積感染者数） $R$ の割合です。これから、基本再生産数 $R_0$ の値によって、最終的にどれくらいの割合で感染者が出るか、どれくらいの収束期間がかかるかを読み取ることができます。ただし、どの図も $R_0$ が一定の場合であることに注意してください。これらの図から、もし $R_0$ の値を順次下げていくことができれば、最終的な感染者数の割合は小さくなることが分かります。

図4 基本再生産数 $R_0$ に依存する感染ダイナミックス

図5 基本再生産数 $R_0=2.5$ を $(1-q)R_0$ に低減したときの効果

$(22)\quad J(R_0)=(I(t_1)/N-i(t_1))^2+\cdots+(I(t_n)/N-i(t_n))^2$

を最小化するようなパラメータ $R_0$ の値 $R_0^*$ を、数値最適化の手法を用いて求めます。

$(23)\quad \underbrace{1.08}_{R_0^*}=1+\lambda_0\times\underbrace{1}_{T}\quad\Rightarrow\quad \lambda_0=0.08}$

となります。これをたとえば感染性時間 $T=1/\gamma=1/0.16=6.25$ を用いて

$(24)\quad R_0=1+\underbrace{0.08}_{\lambda_0}\times\underbrace{6.25}_{T}=1.5}$

のように修正すると、公表値と一致します。

5 おわりに
「なぜ外出自粛は必要か」の問いには、感染爆発を防ぐためには基本再生産数を小さく抑え込むことが必要で、そのための間接的手段として外出自粛やロックダウンが考えられるからと答えられます。感染ダイナミックスの簡易モデル（SIRモデル）は発散から収束に切替る振舞いをうまくモデル化しています。ただ、基本再生産数というたった１つのパラメータで、すべての状況を反映させようとしているところには自ずと限界があります。それでも感染ダイナミクスの本質を捉えていることは間違いないと思います。

「８割削減の妥当性」についてですが、注3で述べたように、基本再生産数 $R_0=2.5$ の場合、図5から「基本再生産数を８割削減することが必要で、２割削減では効果がない」という主張は理解できます。ところが、4/11放映のNHKスペシャルより、基本再生産数 $R_0$ は大体1.5程度であることが表明されています。この場合、図5は違ったものになります。また、5/1の専門家会議では図3のように実効再生産数は4/10の時点で1を切っていることが表明されています。いつの時点で、どのくらいの基本再生産数であるかは、感染対策を議論する上で極めて重要です。もし4/7の時点で基本再生産数が１を切っているのなら、非常事態宣言は医療崩壊を防ぐ意味合いが強かったのかもしれません。

感染流行を制御する決め手は基本再生産数 $R_0$ をどう調整していくかにかかっています。これをSIRモデル(18)の非線形最適制御問題または線形パラメータ変動制御問題として定式化し、その解として望ましい $R_0$ の時系列を求めることが考えられます。もちろん、これが得られたとしても、どのような行動変容に結び付ければよいかは不明です。ただ、厳しい外出自粛をいつまでも続けるわけにはいきません。現在収束過程にあると思われる感染者数データから $R_0$ の時系列を求め、比較検討しておくことは、第２波への対応や出口戦略の検討に役立つと思っております。年内には検討結果をご報告したいと思います。

補遺
5/5に小田垣先生が新型コロナウイルスの蔓延に関する一考察を発表されました。現在指摘されている「PCR検査不足」の主張の妥当性を議論する際のベースを与えることのできるきわめて興味深い論文です。以下では、本稿の枠組みの中でフォローアップさせていただきたいと思います。

当該論文では、感染性人口 $I(t)$ を

● $H(t)$ ：未検査感染性人口（市中感染者数）　⇒　筆者の方で記号をIからHに変更
● $Q(t)$ ：検査後隔離状態にある感染性人口（隔離感染者数）

に分けて、次のSIQRモデルが提案されています。

これをSIRモデル

$(2)\quad\left\{\begin{array}{l} \frac{dS(t)}{dt}=-\beta S(t)I(t) \\ \frac{dI(t)}{dt}=-\gamma I(t)+\beta S(t)I(t)\\ \frac{dR(t)}{dt}=\gamma I(t) \end{array}\right.$

と比較します。それぞれの第 $i$ 式を、SIRi、SIQRiと参照します。

$(26)\quad \beta S(t)I(t)=(1-q)\beta S(t)I(t)+q\beta S(t)I(t)$

このとき、SIRモデルの場合の(5)と同様に

$(27)\quad \frac{dH(t)}{dt}=\underbrace{((1-q)\beta S_0-\gamma-p)}_{\lambda_0(q,p)}H(t), \quad H(0)=H_0$

が成り立ち、初期成長率 $\lambda_0(q,p)$ は

$(28)\quad \lambda_0(q,p)=(\gamma+p)((1-q)\underbrace{\frac{\beta S_0}{\gamma+p}}_{\tilde{R}_0}-1)$

と表すことができます。これをSIRモデルの場合の

$(7)\quad \lambda_0=\gamma(\underbrace{\frac{\beta S_0}{\gamma}}_{R_0}-1)$

(27)の解は

$(29)\quad H(t)=\exp(\lambda_0(q,p)t)H_0\rightarrow0\quad(t\rightarrow\infty)$

と表され、初期値 $H_0$ の10%まで下がる時間は( $\lambda_0(q,p)<0$ に注意して)

$(30)\quad \exp(\lambda_0(q,p)t)=0.1\quad\Rightarrow\quad t=\frac{\ln(0.1)}{|\lambda_0(q,p)|}=\frac{2.3}{|\lambda_0(q,p)|}$

で求められます。

当該論文の重要な貢献は、初期成長率を

$(31)\quad \lambda_0=\beta S_0-\gamma\quad\Rightarrow\quad \lambda_0(q,p)=(1-q)\beta S_0-(\gamma+p)$

これらは次を(31)と(30)に代入して確認できます。

$(32)\quad \beta S_0=0.126$ 、 $\gamma=0.03$ 、 $p=0.096$

たとえば、Case 1)については

当該論文では、基本再生産数を

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.126}{0.03+0.096}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.051}{0.03+0.096}=0.405$

注5：本稿の計算では、基本再生産数は

感染増大期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.128}{0.04}=3.191$
感染減少期間： $R_0=\frac{\beta S_0}{\gamma}=\frac{0.022}{0.04}=0.545$

となります。そこで当該論文と同様の計算をしてみますと

感染頂上期間： $\tilde{R}_0=\frac{\beta S_0}{\gamma+p}=\frac{0.128}{0.04+0.088}=1$
感染減少期間： $\tilde{R}_0=\frac{(1-q)\beta S_0}{\gamma+p}=\frac{0.022}{0.04+0.088}=0.545$

となります。これから削減率 $q=1-0.022/0.128=0.829$ 、隔離率 $p=0.088$ 、 $1/p=11.4$ が分かります。これらを用いると次の計算結果を得ます。

謝辞
最後に、本サイトでまとめた感染シミュレーションは、３月末からわたって検討してきました。その間多くの知人に補遺を除く本サイトを見ていただき、コメントをいただき、その都度見直すことができて、大変感謝申し上げております。本来ならばお名前を挙げて謝意を表すべきところですが、まだ最終原稿となっていない今の段階では、文責を筆者のみに帰すため、遠慮させていただきます。

5/17/2020 梶原宏之　kajiwara@cacsd2.sakura.ne.jp

追記
有料サイトのため、全文を読むことができませんが、次の発表がありました。 $\gamma$ の値は0.04（感染性期間25日）は長すぎるかもしれません。

「コロナ感染は発症前から発症5日が最多、6日以後はほとんど感染しない！」
台湾では2020年5月14日現在、新型コロナウイルス感染症（COVID-19）患者440人、死者7人で感染制御に成功しています。その台湾衛生福利部疾病管制署（台湾CDC）からJAMA Internal Medicine（2020年5月1日オンライン版）に「台湾のCOVID-19感染ダイナミクス」の論文が出ました。100人の感染確定者と2,761人の接触者について追跡（contact tracing）を行い、時期による感染力の違いが明らかになりました。最重要点は次の3つです。接触者追跡は発症4日前までさかのぼれ、COVID-19は発症前から発症5日までに感染リスクがあり、6日以降はない！COVID-19の二次感染までの期間（serial interval）は４～５日、SARSは8.4日。
西伊豆健育会病院病院長　仲田和正
2020年05月14日 16:50

水中線状構造物：閉ループ系

投稿日時: 2018年8月15日投稿者: cacsd

● $\displaystyle{\bf Derivation \#6}$
$\displaystyle{ \left[\begin{array}{cc|cc} M_{qq}+M_{qq}^a & M_{qp} & 0 & 0\\ M_{pq} & M_{pp}+M_{pp}^b & 0 & 0\\\hline -c_0 I_N & 0 & I_N & 0\\ 0 & -d_0 I_N & 0 & I_N \end{array}\right] \left[\begin{array}{c} \ddot{q} \\ \ddot{p} \\\hline \ddot{q}_a \\ \ddot{p}_b \end{array}\right]}$
$\displaystyle{+\left[\begin{array}{cc|cc} D_{qq}+D_{qq}^a & D_{qp} & 0 & 0\\ D_{pq} & D_{pp}+D_{qq}^b & 0 & 0\\\hline 0 & 0 & {\cal D}_{qq}^a & 0 \\ 0 & 0 & 0 & {\cal D}_{pp}^b \end{array}\right] \left[\begin{array}{c} \dot{q} \\ \dot{p} \\\hline \dot{q}_a \\ \dot{p}_b \end{array}\right]}$
$\displaystyle{ +\left[\begin{array}{cc|cc} K_{qq} & 0 & -K_{q_aq_a} & -K_{q_ap_b} \\ 0 & K_{pp} & -K_{p_bq_a} & -K_{p_bp_b} \\\hline 0 & 0 & c_2 I_N & 0\\ 0 & 0 & 0 & d_2 I_N \end{array}\right] \left[ \begin{array}{c} q \\ p \\\hline q_a \\ p_b \end{array} \right] = \left[\begin{array}{c} f_1 \\ 0 \\\hline 0\\ 0 \end{array}\right] }$

$\displaystyle{M_{qq}+M_{qq}^a=[m^y_{ij} + \displaystyle\sum_{l,k=1}^NF^y_{ilkj}q_lq_k]_{i,j=1,\cdots,N}+\chi_2I_N}$
$\displaystyle{M_{qp}=[\displaystyle\sum_{l,k=1}^N N^y_{ilkj} q_l p_k]_{i,j=1,\cdots,N}}$
$\displaystyle{M_{pq}=[\displaystyle\sum_{l,k=1}^N N^z_{ilkj} p_l q_k]_{i,j=1,\cdots,N}}$
$\displaystyle{M_{pp}+M_{qq}^b=[m^z_{ij} + \displaystyle\sum_{l,k=1}^NF^z_{ilkj}p_lp_k]_{i,j=1,\cdots,N}+\chi_2I_N}$

for i=1:N, for j=1:N
Mqq(i,j)=m_y(i,j);
Mqp(i,j)=0;
Mpq(i,j)=0;
Mpp(i,j)=m_z(i,j);
for l=1:N, for k=1:N
Mqq(i,j)=Mqq(i,j)+F_y(i,l,k,j)*q(l)*q(k);
Mqp(i,j)=Mqp(i,j)+N_y(i,l,k,j)*q(l)*p(k);
Mpq(i,j)=Mpq(i,j)+N_z(i,l,k,j)*p(l)*q(k);
Mpp(i,j)=Mpp(i,j)+F_z(i,l,k,j)*p(l)*p(k);
end, end, end, end
MM=eye(4*N,4*N);
MM(1:2*N,1:2*N)=[Mqq+chi2*eye(N,N) Mqp;Mpq Mpp+chi2*eye(N,N)];
MM(2*N+1:3*N,1:N)=-c_0*eye(N,N);
MM(3*N+1:4*N,N+1:2*N)=-d_0*eye(N,N);

$\displaystyle{D_{qq}+D_{qq}^a=[c^y_{ij} +\chi_3\sum_{j=1}^N\int_0^1|v-\dot{\eta}|\phi_i\phi_jd\xi + \displaystyle\sum_{l,k=1}^ND^y_{ilkj}q_lq_k+ \sum_{l,k=1}^NE^y_{ilkj}q_l\dot{q}_k]_{i,j=1,\cdots,N}}$
$\displaystyle{D_{qp}=[\displaystyle\sum_{l,k=1}^NL^y_{ilkj}q_lp_k+\displaystyle\sum_{l,k=1}^NM^y_{ilkj}q_l\dot{p}_k]_{i,j=1,\cdots,N}}$
$\displaystyle{D_{pq}=[\displaystyle\sum_{l,k=1}^NL^z_{ilkj}p_lq_k+\displaystyle\sum_{l,k=1}^NM^z_{ilkj}p_l\dot{q}_k]_{i,j=1,\cdots,N} }$
$\displaystyle{D_{pp}+D_{qq}^b=[c^z_{ij} +\chi_3\sum_{j=1}^N\int_0^1|\dot{\eta}|\psi_i\psi_jd\xi + \displaystyle\sum_{l,k=1}^ND^z_{ilkj}p_lp_k+ \displaystyle\sum_{l,k=1}^NE^z_{ilkj}p_l\dot{p}_k]_{i,j=1,\cdots,N} }$
$\displaystyle{{\cal D}_{qq}^a=[c_1\displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{ak}q_{al}]_{i,j=1,\cdots,N}-c_1I_N}$
$\displaystyle{{\cal D}_{pp}^b=[d_1\displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{bk}p_{bl}]_{i,j=1,\cdots,N}-d_1I_N}$

for i=1:N, for j=1:N
Dqq(i,j)=c_y(i,j)+chi3*c_simp*(abs(V-dy).*f_y(:,i).*f_y(:,j));
Dqp(i,j)=0;
Dpq(i,j)=0;
Dpp(i,j)=c_z(i,j)+chi3*c_simp*(abs(dy).*f_z(:,i).*f_z(:,j));
Dqqa(i,j)=0;
Dppb(i,j)=0;
for l=1:N, for k=1:N
Dqq(i,j)=Dqq(i,j)+D_y(i,l,k,j)*q(l)*q(k)+E_y(i,l,k,j)*q(l)*dq(k);
Dqp(i,j)=Dqp(i,j)+L_y(i,l,k,j)*q(l)*p(k)+M_y(i,l,k,j)*q(l)*dp(k);
Dpq(i,j)=Dpq(i,j)+L_z(i,l,k,j)*p(l)*q(k)+M_z(i,l,k,j)*p(l)*dq(k);
Dpp(i,j)=Dpp(i,j)+D_z(i,l,k,j)*p(l)*p(k)+E_z(i,l,k,j)*p(l)*dp(k);
Dqqa(i,j)=Dqqa(i,j)+c_1*S_y(i,j,k,l)*qa(k)*qa(l);
Dppb(i,j)=Dppb(i,j)+d_1*S_z(i,j,k,l)*pb(k)*pb(l);
end, end, end, end
DD=eye(4*N,4*N);
DD(1:2*N,1:2*N)=[Dqq Dqp;Dpq Dpp];
DD(2*N+1:3*N,2*N+1:3*N)=Dqqa-c_1*eye(N,N);
DD(3*N+1:4*N,3*N+1:4*N)=Dppb-d_1*eye(N,N);

$\displaystyle{K_{qq}=[k^y_{ij} + \displaystyle\sum_{l,k=1}^NB^y_{ilkj}q_kq_l+ \displaystyle\sum_{l,k=1}^NH^y_{ilkj}p_kp_l]_{i,j=1,\cdots,N}}$
$\displaystyle{K_{pp}=[k^z_{ij} + \displaystyle\sum_{l,k=1}^NB^z_{ilkj}p_kp_l+ \displaystyle\sum_{l,k=1}^NH^z_{ilkj}q_kq_l]_{i,j=1,\cdots,N}}$
$\displaystyle{K_{q_aq_a}=[\chi_4 \displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{vk}q_{vl}]_{i,j=1,\cdots,N}}$
$\displaystyle{K_{q_ap_b}=[-\chi_5 \displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{vk}\dot{q}_{l}]_{i,j=1,\cdots,N}}$
$\displaystyle{K_{p_bq_a}=[\chi_4 \displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{vk}p_{vl}]_{i,j=1,\cdots,N}}$
$\displaystyle{K_{p_bp_b}=[\chi_5 \displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{vk}\dot{p}_{l}]_{i,j=1,\cdots,N}}$

for i=1:N, for j=1:N
Kqq(i,j)=k_y(i,j);
Kpp(i,j)=k_z(i,j);
Kqqa(i,j)=0;
Kqpa(i,j)=0;
Kpqb(i,j)=0;
Kppb(i,j)=0;
for l=1:N, for k=1:N
Kqq(i,j)=Kqq(i,j)+B_y(i,l,k,j)*q(l)*q(k)+H_y(i,l,k,j)*p(l)*p(k);
Kpp(i,j)=Kpp(i,j)+B_z(i,l,k,j)*p(l)*p(k)+H_z(i,l,k,j)*q(l)*q(k);
Kqqa(i,j)=Kqqa(i,j)+chi4*S_y(i,j,k,l)*q(k)*q(l);
Kqpa(i,j)=Kqpa(i,j)-chi5*S_y(i,j,k,l)*q(k)*dp(l);
Kpqb(i,j)=Kpqb(i,j)+chi4*S_z(i,j,k,l)*p(k)*p(l);
Kppb(i,j)=Kppb(i,j)+chi5*S_z(i,j,k,l)*p(k)*dp(l);
end, end, end, end
KK=eye(4*N,4*N);
KK(1:N,1:N)=Kqq;
KK(N+1:2*N,N+1:2*N)=Kpp;
KK(1:2*N,2*N+1:4*N)=-[Kqqa Kqpa;Kpqb Kppb];
KK(2*N+1:3*N,2*N+1:3*N)=c_2*eye(N,N);
KK(3*N+1:4*N,3*N+1:4*N)=d_2*eye(N,N);

$\displaystyle{f_1=[\chi_3\int_0^1|v-\dot{\eta}|\phi_iv_id\xi]_{i=1,\cdots,N}}$

f1=zeros(N,1);
for i=1:N
f1(i)=chi3*c_simp*(abs(V-dy).*f_y(:,i).*V);
end

●for $i,j,k,l=1,\cdots,N$

$\displaystyle{m^y_{ij}=\int_0^1\phi_i\phi_jd\xi+ \Gamma_m\phi_i(1)\phi_j(1) }$
$\displaystyle{m^z_{ij}=\int_0^1\psi_i\psi_jd\xi+ \Gamma_m\psi_i(1)\psi_j(1) }$

for i=1:N, for j=1:N
m_y(i,j)=c_simp*(f_y(:,i).*f_y(:,j))+G_m*f_y(n_d+1,i)*df_y(n_d+1,j);
end, end
m_z=m_y;

$\displaystyle{E^y_{ijk\ell}=\int_0^1 \phi_i\phi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m \int_0^1\phi_i\phi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi}$
$\displaystyle{E^z_{ijk\ell}=\int_0^1 \psi_i\psi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m \int_0^1\psi_i\psi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
E_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_y(:,k).*df_y(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_y(:,k).*df_y(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_y(:,k).*df_y(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_y(:,k).*df_y(:,l))));
end, end, end, end
E_z=E_y;

$\displaystyle{F^y_{ijk\ell}=\int_0^1 \phi_i\phi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m \int_0^1\phi_i\phi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi}$
$\displaystyle{F^z_{ijk\ell}=\int_0^1 \psi_i\psi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m \int_0^1\psi_i\psi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
F_y(i,j,k,l)=E_y(i,j,k,l);
end, end, end, end
F_z=F_y;

$\displaystyle{M^y_{ijk\ell}= \int_0^1 \phi_i\phi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m \int_0^1\phi_i\phi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi }$
$\displaystyle{M^z_{ijk\ell}= \int_0^1 \psi_i\psi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m \int_0^1\psi_i\psi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi }$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
M_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_z(:,k).*df_z(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_z(:,k).*df_z(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))));
end, end, end, end
M_z=M_y;

$\displaystyle{N^y_{ijk\ell}= \int_0^1 \phi_i\phi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m \int_0^1\phi_i\phi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi }$
$\displaystyle{+\int_0^1 (\phi_i \phi_j''\psi_k''\psi_\ell''+\phi_i \phi_j''\psi_k'\psi_\ell'''+3\phi_i \phi_j'\psi_k''\psi_\ell'''+\phi_i \phi_j'^T\psi_k'\psi_\ell'''' )d\xi}$
$\displaystyle{N^z_{ijk\ell}= \int_0^1 \psi_i\psi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }$
$\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m \int_0^1\psi_i\psi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi }$
$\displaystyle{+\int_0^1 (\psi_i \psi_j''\phi_k''\phi_\ell''+\psi_i \psi_j''\phi_k'\phi_\ell'''+3\psi_i \psi_j'\phi_k''\phi_\ell'''+\psi_i \psi_j'\phi_k'\phi_\ell'''' )d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
N_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_z(:,k).*df_z(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_z(:,k).*df_z(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+c_simp*(f_y(:,i).*ddf_y(:,j).*ddf_z(:,k).*ddf_z(:,l)…
+f_y(:,i).*ddf_y(:,j).*df_z(:,k).*dddf_z(:,l)…
+3*f_y(:,i).*df_y(:,j) .*ddf_z(:,k).*dddf_z(:,l)…
+f_y(:,i).*df_y(:,j) .*df_y(:,k).*ddddf_y(:,l));
end, end, end, end
N_z=N_y;

$\displaystyle{S^y_{ijkl}=\int_0^1 \phi_i\phi_j\phi_k\phi_l d\xi}$
$\displaystyle{S^z_{ijkl}=\int_0^1 \psi_i\psi_j\psi_k\psi_l d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
S_y(i,j,k,l)=c_simp*(f_y(:,i).*f_y(:,j).*f_y(:,k).*f_y(:,l));
end, end, end, end
S_z=S_y;

$\displaystyle{c^y_{ij}= \int_0^1 2 \sqrt{\beta}{u} \phi_i \phi'_j d\xi }$
$\displaystyle{c^z_{ij}= \int_0^1 2 \sqrt{\beta}{u} \psi_i \psi'_j d\xi }$

for i=1:N, for j=1:N
c_y(i,j)=2*sqrt(beta)*c_simp*(U.*f_y(:,i).*df_y(:,j));
end, end
c_z=c_y;

$\displaystyle{k^y_{ij}= \int_0^1\gamma\phi_i\phi'_jd\xi+\gamma\Gamma_p\phi_i(1)\phi'_j(1) }$
$\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \phi_i \phi_j''d\xi }$ $\displaystyle{+\int_0^1 {u^2} \phi_i\phi''_j d\xi }$ $\displaystyle{+ \int_0^1 \phi_i \phi''''_j d\xi}$
$\displaystyle{k^z_{ij}= \int_0^1\gamma\psi_i\psi'_jd\xi+\gamma\Gamma_p\psi_i(1)\psi'_j(1) }$
$\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \psi_i \psi_j''d\xi }$ $\displaystyle{+\int_0^1 {u^2} \psi_i\psi''_j d\xi }$ $\displaystyle{+ \int_0^1 \psi_i \psi''''_j d\xi}$

for i=1:N, for j=1:N
k_y_0(i,j)=c_simp*(gamma*f_y(:,i).*df_y(:,j))+gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j)…
-c_simp*(gamma*(1-xi+G_p))+c_simp*(f_y(:,i).*ddddf_y(:,j));
end, end

for i=1:N, for j=1:N
k_y(i,j)= k_y_0(i,j)+sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j))…
+c_simp*(U.^2.*f_y(:,i).*ddf_y(:,j));
end, end
k_z=k_y;

$\displaystyle{B^y_{ijk\ell}={1\over 2}\gamma \int_0^1 \phi_i\phi'_j\phi'_k \phi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \phi_i(1)\phi'_j(1) \phi'_k(1) \phi'_\ell(1)}$
$\displaystyle{-{3\over 2} \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \phi_i \phi_j'' \phi'_k \phi'_\ell d\xi }$
$\displaystyle{+\int_0^1 \sqrt{\beta}{\dot{u}}\phi_i \phi_j'' \int_\xi^1 \phi'_k \phi'_\ell d\xi d\xi +\int_0^1 {u^2} ( \phi_i\phi_j' \phi'_k \phi''_\ell- \phi_i\phi''_j\int_\xi^1 \phi'_k \phi''_\ell d\xi) d\xi }$
$\displaystyle{+\int_0^1 (\phi_i \phi_j''\phi_k''\phi_\ell''+4\phi_i \phi_j'\phi_k''\phi_\ell'''+\phi_i \phi_j''\phi_k''\phi_\ell'''' )d\xi}$
$\displaystyle{B^z_{ijk\ell}={1\over 2}\gamma \int_0^1 \psi_i\psi'_j\psi'_k \psi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \psi_i(1)\psi'_j(1) \psi'_k(1) \psi'_\ell(1)}$
$\displaystyle{-{3\over 2} \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \psi_i \psi_j'' \psi'_k \psi'_\ell d\xi }$
$\displaystyle{+\int_0^1 \sqrt{\beta}{\dot{u}}\psi_i \psi_j'' \int_\xi^1 \psi'_k \psi'_\ell d\xi d\xi +\int_0^1 {u^2} ( \psi_i\psi_j' \psi'_k \psi''_\ell- \psi_i\psi''_j\int_\xi^1 \psi'_k \psi''_\ell d\xi) d\xi }$
$\displaystyle{+\int_0^1 (\psi_i \psi_j''\psi_k''\psi_\ell''+4\psi_i \psi_j'\psi_k''\psi_\ell'''+\psi_i \psi_j''\psi_k''\psi_\ell'''' )d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
B_y_0(i,j,k,l)=1/2*gamma*c_simp*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l))…
+1/2*gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j).*df_y(n_d+1,k).*df_y(n_d+1,l)…
-3/2*c_simp*(gamma*(1-xi+G_p))…
+c_simp*(f_y(:,i).*ddf_y(:,j).*ddf_y(:,k).*ddf_y(:,l)…
+4*f_y(:,i).*df_y(:,j) .*ddf_y(:,k).*dddf_y(:,l)…
+f_y(:,i).*ddf_y(:,j).*ddf_y(:,k).*ddddf_y(:,l));
end, end, end, end

for i=1:N, for j=1:N, for k=1:N, for l=1:N
B_y(i,j,k,l)= B_y_0(i,j,k,l)…
+3/2*sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l))…
+sqrt(beta)*c_simp*(dU.*f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*df_y(:,l))))…
+c_simp*(U.^2.*(f_y(:,i).*df_y(:,j).*df_y(:,k).*ddf_y(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*ddf_y(:,l)))));
end, end, end, end
B_z=B_y;

$\displaystyle{H^y_{ijk\ell}={1\over 2}\gamma \int_0^1 \phi_i\phi'_j\psi'_k \psi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \phi_i(1)\phi'_j(1) \psi'_k(1) \psi'_\ell(1)}$
$\displaystyle{-{1\over 2} \int_0^1(\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \phi_i \phi_j'' \psi'_k \psi'_\ell d\xi }$
$\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \phi_i \phi_j' \psi'_k \psi''_\ell d\xi }$
$\displaystyle{-\int_0^1 \sqrt{\beta}{\dot{u}}\phi_i \phi_j'' \int_\xi^1 \psi'_k \psi'_\ell d\xi d\xi +\int_0^1 {u^2} ( \phi_i\phi_j' \psi'_k \psi''_\ell - \phi_i\phi''_j\int_\xi^1 \psi'_k \psi''_\ell d\xi) d\xi }$
$\displaystyle{H^z_{ijk\ell}={1\over 2}\gamma \int_0^1 \psi_i\psi'_j\phi'_k \phi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \psi_i(1)\psi'_j(1) \phi'_k(1) \phi'_\ell(1)}$
$\displaystyle{-{1\over 2} \int_0^1(\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \psi_i \psi_j'' \phi'_k \phi'_\ell d\xi }$
$\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)}$ $\displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) }$ $\displaystyle{) \psi_i \psi_j' \phi'_k \phi''_\ell d\xi }$
$\displaystyle{-\int_0^1 \sqrt{\beta}{\dot{u}}\psi_i \psi_j'' \int_\xi^1 \phi'_k \phi'_\ell d\xi d\xi +\int_0^1 {u^2} ( \psi_i\psi_j' \phi'_k \phi''_\ell - \psi_i\psi''_j\int_\xi^1 \phi'_k \phi''_\ell d\xi) d\xi }$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
H_y_0(i,j,k,l)=1/2*gamma*c_simp*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l))…
+1/2*gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j).*df_y(n_d+1,k).*df_y(n_d+1,l)…
-1/2*c_simp*(gamma*(1-xi+G_p).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l))…
-c_simp*(gamma*(1-xi+G_p).*f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l));
end, end, end, end

for i=1:N, for j=1:N, for k=1:N, for l=1:N
H_y(i,j,k,l)= H_y_0(i,j,k,l)…
+1/2*sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l)…
+sqrt(beta)*(1-xi).*f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l)…
+sqrt(beta)*f_y(:,i).*ddf_y(:,j).*(simp_b*(df_z(:,k).*df_z(:,l))))…
+c_simp*(U.^2.*(f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*ddf_y(:,l)))));
end, end, end, end
H_z=H_y;

$\displaystyle{D^y_{ijk\ell}=\int_0^1 2 \sqrt{\beta}{u} ( \phi_i\phi_j' \phi'_k \phi'_\ell- \phi_i\phi''_j\int_\xi^1 \phi'_k \phi'_\ell d\xi) d\xi}$
$\displaystyle{D^z_{ijk\ell}=\int_0^1 2 \sqrt{\beta}{u} ( \psi_i\psi_j' \psi'_k \psi'_\ell- \psi_i\psi''_j\int_\xi^1 \psi'_k \psi'_\ell d\xi) d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
D_y(i,j,k,l)=2*sqrt(beta)*c_simp*(U.*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*df_y(:,l)))));
end, end, end, end
D_z=D_y;

$\displaystyle{L^y_{ijk\ell}=\int_0^1 2 \sqrt{\beta}{u} ( \phi_i\phi_j' \psi'_k \psi'_\ell- \phi_i\phi''_j\int_\xi^1 \psi'_k \psi'_\ell d\xi) d\xi}$
$\displaystyle{L^z_{ijk\ell}=\int_0^1 2 \sqrt{\beta}{u} ( \psi_i\psi_j' \phi'_k \phi'_\ell- \psi_i\psi''_j\int_\xi^1 \phi'_k \phi'_\ell d\xi) d\xi}$

for i=1:N, for j=1:N, for k=1:N, for l=1:N
L_y(i,j,k,l)=2*sqrt(beta)*c_simp*(U.*(f_y(:,i).*df_y(:,j).*df_z(:,k).*df_z(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_z(:,k).*df_z(:,l)))));
end, end, end, end
L_z=L_y;