感染制御

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の時系列は8割削減と5割削減を適当に設定しています(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を5割削減と2割削減を適当に設定して、向こう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)を1次遅れのダイナミックスを通して間接的に操作できるように政策変数(仮想的な操作変数)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) }

これは2つの変動パラメータ\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} }

いま、2つの変動パラメータを適当な値に\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) }

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

\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.85i^*=0.15T_q=0.2としています。時刻t=3ではすでに人口の15%(0.15N)の感染者がいますので、強力な自粛が必要となり、qの値は1を超える大きな値となっています。そのためシミュレータ上では飽和要素を入れてこれを制限しています。結果的にいきなり10割削減相当、いわゆるロックダウンの自粛となっています。しかし最終的な累積感染者数は9割から3割程度に抑えられています。

図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) }

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

\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) }

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

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

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

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

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モデルの第2式と第3式は、SIRモデルの第2式を、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.

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

図11 正規化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モデル(全国版)に基づく感染予測

図12 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)を1次遅れのダイナミックスを通して間接的に操作できるようにそれぞれの政策変数(仮想的な操作変数)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}

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

(未定稿)