几个常见的数学模型

房室模型

房室模型通过将研究的整体分为各个部分，然后在研究各个部分之间的转移来达成求解模型的目的。 通常，房室之间的转移可由一组微分方程描述。

药代动力学模型

药代动力学主要研究药物的释放、吸收、分布、代谢和排泄问题，后两个阶段也可统称排除。

药代动力学模型中最常用的为二室模型，这种模型中只存在两个房室：中心室和周边室。 中心室表示血液富集的区域，如血管、心、肺等器官； 周边室表示血液较少的区域，如肌肉等。

对这个模型，通常有以下假设：

1. 两个房室的容积，即血液的体积，保持不变；
2. 药物从一个房室转移到其他阶段的速率正比于药物在该室内的浓度；
3. 只有中心室与体外有药物交换，药物通过释放进入中心室，然后从中心室排除，其他浓度变化忽略不计。

从而，中心室有三个参数$c_1(t), x_1(t), V_1$，周边室也有三个参数$c_2(t), x_2(t), V_2$，分别表示药物浓度、药量和容积。 此外，还有一个函数$f(t)$表示给药速率，常数$k_{ij}$表示阶段之间的转移常数。 我们可以列出以下方程： \(\begin{cases} \dot x_1 (t) &= -k_{12}x_1(t) - k_{13} x_1(t) + k_{21} x_2(t) + f_0(t) \\ \dot x_2 (t) &= k_{12} x_1(t) - k_{21} x_2(t) \end{cases}\) 除此之外，还有： \(x_1(t) = c_1(t) V_1, \; x_2(t) = c_2(t) V_2\) 代入原方程，可得： \(\begin{cases} \dot c_1(t) &= - (k_{12} + k_{13}) c_1(t) + \frac{V_2}{V_1} k_{21} c_2(t) + \frac{f_0(t)}{V_1} \\ \dot c_2(t) &= \frac{V_1}{V_2} k_{12} c_1(t) - k_{21} c_2(t) \end{cases}\)

这是简单的非齐次一阶线性方程组，容易求解。

值得注意的是，给药速率$f_0(t)$对这组微分方程的解有非常大的影响，这意味着欲研究不同的给药方式，就需要给出合适的给药速率函数。

传染病模型

我们使用一系列微分方程模型来对传染病的传播进行建模，其中最简单而常用的两种模型为SIS模型和SIR模型。 这些模型都可以看作房室模型。

SI模型

我们首先研究最简单的模型——SI模型，作为剩下两种的铺垫。

SI模型把人口分为两个部分：易感染者（Susceptibles）和已感染者（Infectives）。

该模型中，我们假设总人数$N$不变，而两者的比例分别为$s(t)$和$i(t)$。 除此之外，我们还假设每一个病人每天能与$\lambda$接触并导致他们患病，这一比例称为日接触率。

那么该模型可由以下这组方程表示： \(\left\{ \begin{aligned} \frac{\mathrm d i}{\mathrm d t} &= \lambda i(t) s(t) \\ i(t) + s(t) &= 1 \\ i(0) &= i_0 \end{aligned} \right.\) 使用第二个等式代换第一个等式中的$s(t)$，即可得： \(\frac{\mathrm di}{\mathrm d t} = \lambda i(t) (1-i(t))\) 具有这种方程的模型称为逻辑斯蒂模型，是我们下一部分研究的重点。

这个方程具有解析解，其解如下： \(i(t) = \frac{1}{1 + \left(\frac{1}{i_0}-1\right)e^{-\lambda t}}\)

不难注意到$i(\infty) = 1$，这意味着所有人最终都会被感染。 考虑到这个系统中不存在治愈这种可能，这个结果是符合直观的。

SIS模型

SIS模型在SI模型基础之上建立了从已感染者回到易感染者的通路，这可以用来对感染并治愈后对疾病的免疫力没有显著提高的疾病进行建模。

SIS模型在前一个模型的基础之上，假设病人每天被治愈的比例为$\mu$，其方程可写为： \(\left\{ \begin{aligned} \frac{\mathrm d i}{\mathrm d t} &= \lambda i (1-i) - \mu i = - \lambda i [i - (1 - \frac{1}{\sigma})] \\ i(0)&= i_0 \end{aligned} \right.\) 其中$\sigma = \frac{\lambda}{\mu}$，表示一个感染期内每个病人的有效接触人数，称为接触数。

这个方程的解析解求解比较困难，但定性的分析不难给出： \(i(\infty) = \begin{cases} 1 - \frac{1}{\sigma} &, \sigma > 1 \\ 0 &, \sigma \le 1 \end{cases}\) 这意味着该疾病能否得到控制与其接触数密切相关。

SIR模型

SIR模型在SI模型的基础上新增了移出者（Removed），用来表示被移出感染系统的人数。 这一移出可以是被治愈并获得免疫者，也可能是病死者。

我们维持此前的定义不变，只是将治愈率$\mu$重命名为移出率，其方程可写为： \(\left\{ \begin{aligned} \frac{\mathrm d i}{\mathrm d t} &= \lambda s i - \mu i \\ \frac{\mathrm d s}{\mathrm d t} &= -\lambda s i \\ i(0) &= i_0 \\ s(0) &= s_0 \end{aligned} \right.\)

我们可以使用这个方程来画出其相轨线： \(\frac{\mathrm d i}{\mathrm d s} = \frac{1}{\sigma s} - 1\)

这个模型不具有解析解，但是平衡点和相轨线分析表明，其趋势与$s_0 \cdot \sigma$密切相关。 若$s_0 \cdot \sigma \ge 1$，那么$i(t)$先上升再下降至零； 若$s_0 \cdot \sigma < 1$，那么$i(t)$单调下降至零，这意味者传染病不会蔓延。

因此，$s_0 \sigma$，即初始易感比率与接触数之积，是预测传染病传播的一项重要数据，称为基本再生数。 常见传染性疾病的基本再生数在一至十之间。 据估计，2019-nCoV的基本再生数约在2.3-5之间。

这项数据为预防传染病蔓延提供了一些方法：

1. 降低$s_0$。 注意到易感比例、感染比例和移出比例之和为一，而初始时刻感染比例通常非常小，可视为零，那么通过接种疫苗等方式增加移出比例可降低$s_0$。
2. 降低$\sigma$。 这可通过降低$\lambda$或提高$\mu$来实现。

逻辑斯蒂模型

逻辑斯蒂模型是一大类解具有逻辑斯蒂函数的形式的模型的总成。 逻辑斯蒂函数，是一种S型函数，可用来表示一个增长既收到其自己的促进，又收到自己的抑制的物理量随时间的变化，如生物种群等。

逻辑斯蒂模型的微分方程通常具有以下形式： \(\frac{\mathrm d x}{\mathrm d t} = r x (1 - \frac{x}{x_m})\) 其通解为： \(x(t) = \frac{x_m}{1 + (\frac{x_m}{x_0} - 1) e^{-rt}}\)

Leslie模型

逻辑斯蒂模型可较好地描述种群的增长，但是却无法给出更多人口性别或年龄分层的数据。 为了更好地确定人口的年龄分组，我们使用Leslie模型。 这个模型是一个差分模型，并非逻辑斯蒂模型。

假设我们将种群分为$n$个年龄层，其中年龄最小的人数为$x_1$。 Leslie模型使用几个差分方程来估计种群数量： \(x_1(k+1) = \sum_{i=1}^n b_i x_i(k), \; x_{i+1}(k+1) = s_i x_i(k)\) 第一个方程表示各个年龄分组产下的新个体，全部进入年龄最小的年龄组内，$b_i$称为出生率。 第二个方程表示年龄分组的成长，$s_i$称为存活率。

若我们把种群数写成向量的形式： \(x(k) = \begin{bmatrix} x_1(k) & \cdots & x_n(k) \end{bmatrix}^\top\) 那么这个差分方程可由一个矩阵与其相乘表示： \(L = \begin{bmatrix} b_1 & b_2 & \cdots & b_{n-1} & b_n \\ s_1 & 0 & \cdots & 0 & 0 \\ 0 & s_2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & s_{n-1} & 0 \end{bmatrix}, x(k+1) = L x(k)\) 矩阵$L$称为Leslie矩阵。

Leslie模型的稳态分析

特征值是对矩阵进行分析的有力手段，因此我们也借助它对Leslie矩阵进行分析。 设其最大单特征值为$\lambda_1$，通过简单的计算可知其特征向量可写为： \(x^* = \begin{bmatrix} 1, \frac{s_1}{\lambda_1}, \frac{s_1 s_2}{\lambda_1^2}, \cdots, \frac{s_1 s_2 \cdots s_{n-1}}{\lambda_1^{n-1}} \end{bmatrix}^\top\)

不妨设该矩阵可对角化： \(L = P\times \mathrm{diag}(\lambda_1, \cdots, \lambda_n) \times P^{-1}\) 注意到所有其他特征值均小于$\lambda_1$，我们有： \(\lim_{k \to \infty} \frac{x(k)}{\lambda_1^k} = P \times \mathrm{diag}(1, 0, \cdots, 0) \times P^{-1} \times x(0)\) 又因为$P$的第一列正是第一个特征向量$x^*$，我们可以认为： \(x(\infty) = c x^*\) 其中$c$是一个由Leslie矩阵和$x(0)$决定的常量。

注意到特征向量$x^*$与初态无关，这个分布向量称为稳定分布。 设$\lambda = \lambda_1$，当$k$充分大时，可以认为： \(x(k) = c \lambda^k x^*\)

此外， \(x(k+1) = \lambda x(k) \implies x_i(k+1) = \lambda x_i(k)\) 因此$\lambda$也被称为固有增长率。

相互影响的种群模型

利用逻辑斯蒂模型，可对相互影响的两个种群的数量进行建模。

竞争的种群

设甲乙两个种群相互竞争，则两者会争夺同一资源，从而互相对对方的种群增长产生阻滞。

设甲乙两个种群的数量服从逻辑斯蒂模型，分别为$x_1, x_2$，其固有增长率为$r_1, r_2$，容量为$N_1, N_2$，那么两者独立时，甲的微分方程可写为： \(\dot x_1(t) = r_1 x_1 \left( 1 - \frac{x_1}{N_1} \right)\) 乙的与其类似。

现在，若考虑两者之间的竞争。 对甲而言，显然，其增长速度收到乙种群数量的限制。 我们设其限制正比于乙的种群数量，那么微分方程可写为： \(\dot x_1(t) = r_1 x_1 \left(1 - \frac{x_1}{N_1} - \sigma_1 \frac{x_2}{N_2} \right)\) 这里的$\sigma_1$可以解释为，单位数量的乙（相对$N_2$而言）消耗的资源数量是单位数量甲（相对$N_1$而言）的$\sigma_1$倍。 同理，乙的微分方程可写为： \(\dot x_2(t) = r_2 x_2 \left(1 - \frac{x_2}{N_2} - \sigma_2 \frac{x_1}{N_1} \right)\)

若$\sigma_1, \sigma_2$描述的是两个种群争夺同一资源的能力，那么通常而言二者的乘积应当为一。 但是，这种情况比较理想，因此我们希望研究两个变量独立的情况。

结合平衡点的局部稳定性和相轨线分析，系统存在四种平衡状态，其中三种全局稳定：

1. $P = (N_1, 0)$，$\sigma_1 < 1, \sigma_2 > 1$时稳定。
2. $P = (0, N_2)$，$\sigma_1 > 1, \sigma_2 < 1$时稳定。
3. $P = (c_1N_1, c_2N_2)$，$\sigma_1, \sigma_2 < 1$时稳定。
4. $P = (0, 0)$，不稳定。

依存的种群

现在，设种群甲可以独立存在，种群乙和种群甲相互依存，但必须依赖于种群甲才能存在。 那么此时种群甲的微分方程可写为： \(\dot x_1(t) = r_1 x_1 \left(1 - \frac{x_1}{N_1} + \sigma_1 \frac{x_2}{N_2} \right)\)

对种群乙，首先，当种群甲不存在时，其自然消亡，有： \(\dot x_2(t) = -r_2 x_2\) 而种群甲的存在对其起促进作用，从而： \(\dot x_2(t) = r_2 x_2 (-1 + \sigma_2 \frac{x_1}{N_1})\) 然而，种群乙对自身的增长也有抑制作用，从而其完整的微分方程为： \(\dot x_2(t) = r_2 x_2 (-1 - \frac{x_2}{N_2} + \sigma_2 \frac{x_1}{N_1})\)

这个模型有三个平衡点，其中两个稳定： \(\begin{array}{ll} P_1 = (N_1, 0)&, \sigma_2 < 1, \sigma_1 \sigma_2 < 1 \\ P_2 = (N_1 \frac{1-\sigma_1}{1 - \sigma_1\sigma_2}, N_2 \frac{\sigma_2-1}{1-\sigma_1\sigma_2})&, \sigma_1 < 1, \sigma_2 > 1, \sigma_1\sigma_2 < 1 \\ P_3 = (0, 0)&, \text{不稳定} \end{array}\)

优化模型

本节主要研究静态优化问题，即最优解是一个数而非函数的优化问题。 优化问题的关键在于设计一个性质良好的目标函数。 这种优化问题通常使用微分法求解。

单变量优化模型：仓储问题

首先我们来解决最简单的一种优化模型：单变量优化模型。

我们将要求解管理学中最简单的一种仓储模型，即经济订货量（Economic Ordering Quantity，EOQ）模型。 假设一座仓库中存有一些日需求量为常数$r$的货品，该产品每次生产消耗$c_1$元，每$T$天生产一次，每次生产$Q$件，每日每件产品的储存费为$c_2$元。 现在假设$r, c_1, c_2$已知，在不发生缺货的情况下，试求$T,Q$使总花费最小。

我们假设变量都是连续的，以此简化求解。 由于不允许缺货，一周期$T$内货物的总需求量必须和生产量相等，因此$Q = rT$。 从而一周期内的储存费用为： \(c_{2,\text{total}} = c_2 \int_0^T q(t) \mathrm d t = c_2 \frac{QT}{2}\) 从而总费用为： \(\tilde C(T) = c_1 + c_2 \frac{QT}{2} = c_1 + c_2 \frac{rT^2}{2}\)

这个函数非常初等，直接求导并使导数等于零即可求出最值。 从而容易知道： \(T = \sqrt{\frac{2c_1}{rc_2}}\) 因此单次最优订货量为： \(Q = \sqrt{\frac{2c_1r}{c_2}}\)

线性规划模型

线性规划模型是数学规划模型中的一个类别。

数学规划模型通常可以形式化地写为： \(\begin{aligned} \min z = f(x) & \quad x = (x_1, \cdots, x_n)^\top\\ s.t. \quad g_i(x) \le 0 & \end{aligned}\) 其中$z=f(x)$称为目标函数，$x$称为决策变量，$g_i(x)$称为约束条件。

若目标函数和所有约束条件都是线性函数，那么这个规划模型就是线性规划模型。

二维的线性规划模型通常可以通过画图完成，然而更高维的问题需要更通用的解决办法。 在此之前，我们首先需要将线性规划模型转化为标准型，以便于求解。

线性规划模型的标准型

线性规划模型通常具有以下形式： \(\begin{aligned} \max (\min) z = c_1 x_1 + c_2 x_2 + \cdots + c_n x_n \\ a_{11} x_1 + \cdots + a_{1n} x_n \; \mathbf{op} \; b_1 \\ a_{21} x_1 + \cdots + a_{2n} x_n \; \mathbf{op} \; b_2 \\ \cdots \qquad \\ a_{m1} x_1 + \cdots + a_{mn} x_n \; \mathbf{op} \; b_m \\ x_j \ge 0 \quad j \in \{ 1, \dots, n \} \end{aligned}\) 其中$\mathbf{op}$表示大于等于、等于或小于等于之一。

而线性规划的标准型具有以下形式： \(\begin{aligned} \min z = c_1 x_1 + c_2 x_2 + \cdots + c_n x_n \\ a_{11} x_1 + \cdots + a_{1n} x_n = b_1 \\ a_{21} x_1 + \cdots + a_{2n} x_n = b_2 \\ \cdots \qquad \\ a_{m1} x_1 + \cdots + a_{mn} x_n = b_m \\ x_j \ge 0 \quad j = 1, \dots, n \end{aligned}\) 注意到在标准型中我们要求只进行最小化、所有约束都是等于约束且所有变量大于等于零。 这让我们能很容易地将线性规划模型写成矩阵的形式： \(\begin{aligned} \min z = CX & \\ s.t. \quad AX = b& \\ X \ge 0 \end{aligned}\)

通常我们使用以下几步将模型化为标准型：

1. 将极大值转化为求极小值：只需要添加负号即可。
2. 将不等式转化为等式：添加松弛变量或剩余变量。 注意到： \(\begin{aligned}x \le b \iff x + x_{\text{rel}} = b \\ x \ge b \iff x - x_{\text{res}} = b\end{aligned} \\ x_\text{rel}, x_{\text{res}} \ge 0\) 通过添加一个非负变量即可将不等式转化为等式。
3. 无约束变量转为约束变量。 注意到： \(x \text{ unconstrained} \iff x^\prime - x^{\prime\prime} \ge 0\) 我们可以使用两个非负变量来替代一个无约束的变量。
4. 负约束转化为非负约束：为变量添加负号即可。

单纯形法

单纯形法使用迭代的方法来求解线性规划问题。

首先，我们定义基本变量和非基本变量。 假设需要求解的线性规划问题的目标函数中具有$n$个变量，而约束中总共出现了$m$个变量。 我们总是有$m \ge n$，多出来的这些变量通常是在转化为标准型时添加的。 那么，这个问题中就有$m-n$个基本变量，$n$个非基本变量。 在求解问题时，通过将一个变量移至等式一侧并使其系数为零，该问题的约束总是可以等价地写为： \(\begin{aligned} x_i = b_i + c_1 x_1 + \cdots + c_m x_m \\ x_j = b_j + c_1 x_1 + \cdots + c_m x_m \\ \dots \dots \qquad \end{aligned}\) 有一些变量只会出现在左侧，而剩下的变量只会出现在右侧。 出现在左侧的变量称为基本变量，出现在右侧的变量称为非基本变量。 通过恒等变形，我们可以让目标函数只通过非基本变量表出。

单纯形算法通过选择并变换这些基本变量来求解线性规划问题，因此在开始时总是需要选择一组基本变量。 如果一个问题的约束只有不等式，那么通过将所有松弛和剩余变量移动到等式的另一侧即可非常方便地构造出一组基本变量。

单纯形算法选择并变换基本变量的操作称为转动（Pivoting）。 这一操作非常类似于高斯消元中的选主元（Pivoting）。 转动操作由一个替出变量和替入变量组成。 顾名思义，该操作就是从非基本变量中选择一个替入变量来替换基本变量中的替出变量。 这一操作非常简单：选择替出变量所在的等式，然后将替入变量移到右边，替出变量移到左边，然后两边同时除以替入变量的系数使其归一，最后重写所有等式，用替入变量替代替出变量即可。 这一系列操作都可以用矩阵的基本行变换表示。

单纯型法的算法步骤如下。

1. 将问题转化为标准型，选择一组初始的基本变量。
2. 判断当前解是否为最优解：计算检验数。 检验数即目标函数中每个变量的系数，若所有系数均大于零，那么这组解就是最优解。 注意目标函数必须表示成非基本变量的线性组合。
3. 选择替出变量和替入变量，进行转动并回到第二步。 检验数最小（即绝对值最大）的变量即为替入变量，替出变量按比值最小的原则进行选择。

我们使用以下例子来介绍这一方法。

我们尝试求解优化问题： \(\begin{aligned} \min z = -5x_1 -2x_2 \\ 30 x_1 + 20 x_2 + x_3 = 160 \\ 5 x_1 + x_2 + x_4 = 15 \\ x_1 + x_5 = 4 \\ x_i \ge 0 \end{aligned}\) 这个问题已经是标准型了，而且变量$x_3, x_4, x_5$的系数都是一，因此我们选择它们作为初始基本变量。

我们使用单纯形表来求解这一问题。 单纯形表的每一列表示一个变量，每一行表示一个基本变量以及其代表的等式。 这个问题的单纯形表为：

	x1	x2	x3	x4	x5	b
x3	30	20	1	0	0	160
x4	5	1	0	1	0	15
x5	1	0	0	0	1	4
-z	-5	-2	0	0	0	0

注意$z$的符号由问题决定，若为最小化问题，则有负号，否则没有。

其检验数，即$z$一行的系数，均为负，因此必须继续迭代。 我们选择检验数最小的变量作为替入变量，即$x_1$。 现在，用这一列的所有非负数除以该行对应的$b$，分别可得：$\frac{160}{3}, 3, 4$。 其中最小者为$x_4$，因此选择它作为替出变量，表格变为：

	x1	x2	x3	x4	x5	b
x3	30	20	1	0	0	160
x1	5	1	0	1	0	15
x5	1	0	0	0	1	4
-z	-5	-2	0	0	0	0

接下来进行基本行变换，首先使$x_1$一列，$x_1$一行的数为一，然后使$x_1$一列的所有数，除了$x_1$一行之外全部变为零。 变换之后的表为：

	x1	x2	x3	x4	x5	b
x3	0	14	0	-6	0	70
x1	1	1/5	0	1/5	0	3
x5	0	-1/5	0	-1/5	1	1
-z	0	-1	0	1	0	15

然后重复迭代，选择$x_2$作为替入变量，忽略系数非正的数并计算比值，选择$x_3$作为替出变量，进行基本行变换，可得：

	x1	x2	x3	x4	x5	b
x2	0	1	1/14	-1/14	0	5
x1	1	0	…	…	0	2
x5	0	0	…	…	1	2
-z	0	0	1/14	2/14	0	20

现在所有检验数非负，因此即为最优解。 最优解在 \(x_1 = 2, x_2 = 5, x_5 = 2\) 时取得，值为$z = -20$。

如果有多个最小的检验数，那么任选一个即可。 如果选择的替入变量一列的系数全部为负，那么该问题无解。 如果出现检验数为零，那么存在无穷多组解。