感染制御 | 制御系CAD

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モデル（線形パラメータ変動モデル）とみなすことができます。