引言与案例
案例1：使用单盘天平称量三个物体
这个经典问题来自Pázman在1986年提出的称量问题。我们有三个物体，它们的真实质量分别为 m_1, m_2, m_3，需要使用一个单盘天平来估计这些质量。
问题的数学模型
当我们用天平称量质量为 m 的物体时，观测值 y 并不是精确的质量值，而是受到两个因素的影响：
这个观测模型包含了实际测量中的两类误差。m_0 代表天平本身存在的系统性偏差，这是一个固定但未知的常数，每次称量都会加上这个偏差。\varepsilon 则是随机测量误差，它在每次测量时都是随机变化的。测量误差满足两个统计性质：期望值为零 E\{\varepsilon\} = 0，这意味着误差没有倾向性；方差为 \text{var}\{\varepsilon\} = \sigma^2，这刻画了误差的波动幅度。这里有一个关键假设：不同次称量的误差是独立同分布的，即每次测量的误差互不影响，且都遵循相同的分布规律。
方法1：单独称量后减去空载
这个方法的核心思想是通过差分消除系统偏差。首先，我们对每个物体进行单独称量，得到三次观测：
每次观测都包含了物体真实质量、系统偏差和随机误差。为了消除 m_0 的影响，我们额外进行一次空载称量（天平上不放任何物体）：
这次空载称量直接测出了系统偏差加上随机误差。通过相减操作，我们构造质量估计量：
系统偏差 m_0 在减法中被完全消除了。现在我们分析这个估计量的统计性质。首先看无偏性：
由于两个误差项的期望都是零，估计量的期望恰好等于真实质量，说明这是一个无偏估计。但是，估计量的方差变大了：
这是因为两次独立测量的误差方差会相加，导致估计的不确定性是单次测量的两倍。更严重的问题在于，不同质量估计之间存在相关性：
这个正相关来源于所有估计量都共同使用了空载测量 y(0)，其中的误差 \varepsilon_0 会同时影响所有三个估计值。这种相关性意味着如果一个质量被高估，其他质量也倾向于被高估。
方法2：满载称量的巧妙设计
方法2同样使用四次称量，但设计更加巧妙。前三次称量与方法1相同：
但第四次称量不是空载，而是将所有三个物体同时放在天平上进行满载称量：
这次测量得到了总质量加上系统偏差和误差。现在，质量估计量的构造方式变得更复杂。对于物体 i，我们使用一个线性组合：
这里 j 和 k 是除了 i 之外的另外两个物体的索引。将观测值代入后展开计算，y(i) + y(0) - y(j) - y(k) 实际上等于 2m_i + m_0 + \varepsilon_i + m_1 + m_2 + m_3 + m_0 + \varepsilon_0 - m_j - m_0 - \varepsilon_j - m_k - m_0 - \varepsilon_k = 2m_i + \varepsilon_i + \varepsilon_0 - \varepsilon_j - \varepsilon_k，因此除以2后系统偏差完全消去，得到无偏估计：
但方差的表现大大改善了：
这个方差只有方法1的一半，因为虽然涉及四次测量的误差，但它们在线性组合中被巧妙地抵消了一部分。最关键的改进是协方差：
不同质量估计之间完全不相关。这是因为每个估计量使用的测量组合都不同，它们之间的误差项通过精心设计实现了统计独立。
因此得出结论：方法2更精确，估计之间没有相关性。这个例子展示了实验设计的核心理念：通过合理安排测量方案，在相同的测量次数下可以获得更好的估计性能，既降低了方差又消除了相关性。
案例2：使用双盘天平称量八个物体
问题背景
现在我们面对一个更复杂的称量问题。有8个物体需要称量，它们的真实质量为 m_i, i = 1, \ldots, 8。与案例1不同，这里使用的是双盘天平，这种天平的工作原理是比较两边的质量差。当我们在左盘放物体1，右盘放物体2时，天平的读数为：
测量误差 \varepsilon_i 服从独立同分布的正态分布 \mathcal{N}(0, \sigma^2)。双盘天平的优势在于它直接测量质量差，不存在系统偏差项 m_0，这是因为两边的系统偏差会相互抵消。
方法1：单独称量并多次平均
最直接的方法是对每个物体单独进行称量。由于双盘天平需要比较，我们可以在一边放物体，另一边空着（或者理解为与零质量比较）。第 i 个物体的单次称量结果为：
基于单次称量的质量估计量为 \hat{m}_i = y(i)，它服从正态分布 \mathcal{N}(m_i, \sigma^2)。这个估计量是无偏的，但方差为 \sigma^2。
为了提高精度，一个自然的想法是对每个物体重复称量8次，然后取平均值。如果对物体 i 进行8次独立称量，得到 y_1(i), y_2(i), \ldots, y_8(i)，那么平均估计量为：
由于这8次测量的误差是独立的，平均后的方差变为原来的 1/8，因此 \hat{m}_i \sim \mathcal{N}(m_i, \sigma^2/8)。但这个方法的代价是需要进行 8 \times 8 = 64 次称量。
方法2：正交设计矩阵
方法2展示了实验设计的精妙之处。我们仍然只进行8次称量，但每次称量时将多个物体巧妙地分配到天平的两边。具体的称量方案如下：
每个方程中，正号表示该物体放在天平的左边，负号表示放在右边。这个设计的核心在于：每个物体在8次称量中，有4次出现在左边（系数为+1），4次出现在右边（系数为-1）。现在看如何从这8次观测中估计 m_1。将前四个观测值加起来，再减去后四个观测值，然后除以8：
将每个 y(i) 的表达式代入，会发现除了 m_1 之外，其他所有质量项的系数在加减后都恰好为零。这是因为对于任意 m_j (j \neq 1)，它在前4个方程中出现的次数（带符号）与在后4个方程中出现的次数正好抵消。最终得到：
由于8个误差项都是独立的，均值为零，方差为 \sigma^2，因此估计量的方差为 \frac{8\sigma^2}{64} = \frac{\sigma^2}{8}，得到 \hat{m}_i \sim \mathcal{N}(m_i, \sigma^2/8)。
这个设计的优越性在于：仅用8次称量就达到了方法1需要64次称量才能达到的精度。方法2背后的数学原理是Hadamard矩阵的正交性，这种设计保证了每个质量的估计都利用了所有8次测量的信息，而且不同质量的估计之间互不干扰。
案例3：药物动力学的双室模型
模型描述
这个案例来自D'Argenio在1981年提出的药物动力学研究。双室模型是描述药物在体内分布和代谢的经典模型。模型包含两个隔室：中心室C（Central compartment）代表血液循环系统，外周室P（Peripheral compartment）代表组织器官。药物通过输入 u(t) 进入中心室，然后在两个隔室之间进行交换，同时从中心室被消除。
系统用两个状态变量描述：x_C(t) 表示中心室的药物量，x_P(t) 表示外周室的药物量。三个速率常数刻画了药物的转运过程：K_{CP} 是从中心室到外周室的转运速率，K_{PC} 是从外周室回到中心室的转运速率，K_{EL} 是从中心室的消除速率。观测方程描述我们实际能测量的量，即中心室的药物浓度：
这里 V 是中心室的分布容积。系统的动力学由以下线性微分方程组描述：
第一个方程表明中心室药物量的变化率等于输入 u(t) 加上从外周室回流的量 K_{PC}x_P(t)，减去转运到外周室的量 K_{CP}x_C(t) 和被消除的量 K_{EL}x_C(t)。第二个方程描述外周室药物量的变化，它只受到两个隔室之间的转运影响。
实际测量时，观测值包含测量误差：
其中 \varepsilon(t) 是独立同分布的测量误差，服从正态分布 \mathcal{N}(0, \sigma^2)，标准差为 \sigma = 0.2 \mu g/ml。
实验设计
实验的目标是估计四个未知参数构成的向量 \mathbf{p} = (K_{CP}, K_{PC}, K_{EL}, V)。输入方案设计为两阶段：在前1分钟内以 u(t) = 75 mg/mn的速率输入药物，然后在1分钟到2分钟期间以较低的 u(t) = 1.45 mg/mn速率继续输入。这种两阶段输入设计能够激发系统的动态响应，帮助识别不同的参数。
为了验证实验设计和参数估计方法，进行了仿真实验，使用的真实参数值为：
这些参数值代表了真实的生理条件，后续的实验设计和参数估计都将围绕如何最有效地从含噪声的浓度测量数据中恢复这些参数展开。
实验设计方案的比较
在药物动力学实验中，实验变量 \xi_i 指的是测量时间点 t_i。实验设计的核心问题是：在整个实验时间范围 1 \leq t_i \leq 720 分钟内，应该在哪些时间点进行浓度测量才能最好地估计参数。
常规方案
常规方案采用均匀间隔的时间点设计：
这8个时间点的选择遵循一个经验规则：在实验初期采样较密（5分钟和10分钟），然后逐渐拉大间隔。这种设计的逻辑是药物浓度在初期变化较快，因此需要更密集的采样来捕捉动态过程，而后期浓度变化缓慢，可以降低采样频率。
最优方案
通过优化算法得到的针对真实参数 \mathbf{p}^* 的最优测量时间方案为：
最优方案呈现出一个显著特征：某些时间点重复出现。具体来说，在1分钟、10分钟、74分钟和720分钟这四个时刻，各进行两次独立测量。这里的关键假设是可以在同一时刻进行独立的测量，即这两次测量的误差是相互独立的。这种设计反映了不同时刻对参数估计的贡献不同：某些关键时刻的信息量特别大，值得投入更多的测量资源。1分钟正好对应输入模式的切换点，74分钟可能对应系统动态响应的特征时刻，而720分钟的终点测量帮助确定长期行为。
蒙特卡罗仿真验证
为了比较两种方案的性能，进行了大规模的蒙特卡罗仿真研究。仿真过程分为三步：首先进行400次独立的随机仿真，每次仿真模拟整个药物动力学过程并在指定时间点生成含噪声的观测数据。对于每种测量方案（常规和最优），都生成400批数据，每批包含8次观测。然后对每批数据使用最小二乘法进行参数估计，得到400组参数估计值。这400组估计值的分布反映了参数估计的统计性质。
估计精度的对比分析
通过仿真得到的参数估计边缘密度分布展示了两种方案的性能差异。图中用实线表示最优方案的密度分布，虚线表示常规方案的密度分布。四个参数的真实值分别标注在图上：K_{CP}^* = 0.066，K_{PC}^* = 0.038，K_{EL}^* = 0.0242，V^* = 30。
从密度曲线可以看出，最优方案的密度分布在所有四个参数上都表现出更尖锐的峰值和更窄的分布范围。对于 \hat{K}_{CP}，最优方案的密度峰高达约18，而常规方案仅约11，且最优方案的分布更集中在真实值附近。对于 \hat{K}_{PC}，两种方案的密度峰值相差更大，最优方案约为24，常规方案约为17。参数 \hat{K}_{EL} 的差异最为显著，最优方案的密度峰值接近170，而常规方案仅约65，说明最优方案对这个参数的估计精度提升非常明显。分布容积 \hat{V} 的估计也显示出类似的改进，最优方案的密度峰值约为0.21，常规方案约为0.04。
密度曲线的宽度反映了估计的方差，更窄的分布意味着更小的方差。最优方案在所有参数上都实现了更小的估计方差，这意味着使用相同数量的测量次数，最优设计能够提供更精确、更可靠的参数估计。这个案例清晰地展示了实验设计优化的价值：不是简单地增加测量次数，而是通过优化测量时间点的选择，就能显著提高参数估计的质量。
实验设计准则
面向估计的实验设计框架
我们考虑基于最大似然估计的面向估计的实验设计问题。在参数估计理论中，最大似然估计量具有良好的渐近性质。当观测次数 N 趋于无穷时，参数估计量 \hat{\mathbf{p}}_{MV}^N 的分布满足渐近正态性：
这个结果说明，当样本量足够大时，估计量的分布近似服从多元正态分布，其均值为真实参数 \mathbf{p}^*（无偏性），协方差矩阵为Fisher信息矩阵的逆 \mathbf{M}^{-1}(\mathbf{p}^*, \boldsymbol{\xi}_1^N)。Fisher信息矩阵 \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}_1^N) 定义为对数似然函数的二阶偏导数的期望：
这个矩阵衡量了观测数据中包含的关于参数的信息量。对数似然函数的梯度 \frac{\partial \log \pi}{\partial \mathbf{p}} 称为得分函数，它描述了观测值对参数的敏感程度。Fisher信息矩阵实际上是得分函数的外积的期望，除以 N 后得到平均每次观测的信息量。
优化准则的选择
实验设计的目标是通过选择合适的实验条件 \boldsymbol{\xi}_1^N 来最小化估计的不确定性。这里的实验变量 \xi_i 可以是多种实验条件，比如观测时间、压力、温度、输入形状等。所有可能的实验条件构成实验域 \mathcal{X}，即 \xi_i \in \mathcal{X}。在很多实际问题中，实验域是有限的离散集合 \mathcal{X} = \{\xi^1, \ldots, \xi^M\}，其中 M 是可选实验条件的数量（通常 M 很大）。优化的准则函数是Fisher信息矩阵的标量函数，它将一个矩阵映射为一个数值，表示实验设计的整体质量。
Kiefer在1974年提出了一个统一的优化准则框架：
这个准则依赖于两个设计参数：幂次 k 和权重矩阵 \mathbf{Q} \succeq \mathbf{O}（半正定）。权重矩阵允许我们对不同参数或参数组合赋予不同的重要性。当 \det \mathbf{M} = 0 时，信息矩阵是奇异的，意味着某些参数组合无法被唯一识别，此时准则值设为无穷大表示该设计不可接受。
A-最优性准则
当选择 k = 1 和 \mathbf{Q} = \mathbf{I}_p（单位矩阵）时，得到A-最优性准则：
矩阵的迹等于其所有特征值之和，对于协方差矩阵 \mathbf{M}^{-1}，迹等于所有参数估计方差的总和。A-最优性的几何意义是最小化渐近置信椭球的轴长之和。在 p 维参数空间中，置信椭球的主轴长度由协方差矩阵的特征值决定，最小化迹相当于让所有方向上的不确定性总和最小。
D-最优性准则
当选择 k = \infty 和 \mathbf{Q} = \mathbf{I}_p 时，当 k \to \infty 时准则函数的极限形式给出D-最优性准则：
由于我们要最小化 \Phi，而D-最优性定义为最大化行列式，所以实际上是最大化 \det \mathbf{M} 或等价地最大化 \log \det \mathbf{M}（因为对数是单调函数）。D-最优性的几何意义是最大化Fisher信息矩阵的特征值乘积，或者等价地最小化渐近置信椭球的体积。椭球的体积正比于 1/\sqrt{\det \mathbf{M}}，因此最大化 \det \mathbf{M} 等价于最小化置信区域的体积，使得参数估计在所有方向上都尽可能精确。
D_s-最优性：部分参数优化
在很多实际问题中，参数向量 \mathbf{p} 可以分为两部分：n_s 个感兴趣的参数和 n_p - n_s 个干扰参数。将参数向量写成 \mathbf{p}^{\mathrm{T}} = (\mathbf{p}_1^{\mathrm{T}}, \mathbf{p}_2^{\mathrm{T}})，其中 \mathbf{p}_1 是 s 维的感兴趣参数向量，\mathbf{p}_2 是干扰参数。相应地，Fisher信息矩阵可以分块表示为：
这里 \mathbf{M}_{11} 是 s \times s 的子矩阵，对应感兴趣参数之间的信息；\mathbf{M}_{22} 是 (n_p - s) \times (n_p - s) 的子矩阵，对应干扰参数之间的信息；\mathbf{M}_{12} 和 \mathbf{M}_{21} 是交叉项。
当我们只关心 \mathbf{p}_1 的估计精度时，需要考虑 \mathbf{M}^{-1} 的相应子块。利用分块矩阵求逆公式，\mathbf{M}^{-1} 可以表示为：
其中各子块的计算公式为：
\mathbf{A}_{11} 给出了在存在干扰参数 \mathbf{p}_2 的情况下，感兴趣参数 \mathbf{p}_1 的渐近协方差矩阵。表达式 \mathbf{M}_{11} - \mathbf{M}_{12}\mathbf{M}_{22}^{-1}\mathbf{M}_{21} 称为Schur补，它表示在消除干扰参数影响后，关于 \mathbf{p}_1 的有效信息矩阵。
D_s-最优性准则的目标是最小化 \mathbf{p}_1 的置信椭球体积，准则函数定义为：
最大化这个行列式等价于最小化 \mathbf{p}_1 的置信区域体积，同时自然地考虑了干扰参数的存在。这个准则在我们只想精确估计部分参数时特别有用，它允许我们将实验资源集中在提高感兴趣参数的估计精度上。
D-最优性的重要性质
D-最优实验设计具有重参数化不变性，这是一个非常重要的性质。假设我们对参数进行变换，从原参数 \mathbf{p} 变换到新参数 \mathbf{p}'，两者通过可微的一一对应关系联系。在新参数下计算的Fisher信息矩阵与原参数下的关系为：
这个等式表明，新旧参数下的行列式相差一个Jacobian矩阵行列式的平方。由于Jacobian行列式 \frac{\partial \mathbf{p}'}{\partial \mathbf{p}^{\mathrm{T}}} 只依赖于参数变换本身，与实验设计 \boldsymbol{\xi} 无关，因此最大化 \det \mathbf{M}(\mathbf{p}', \boldsymbol{\xi}) 等价于最大化 \det \mathbf{M}(\mathbf{p}, \boldsymbol{\xi})。这意味着D-最优设计不依赖于我们选择哪种参数表示方式，无论使用原始参数还是经过变换的参数，得到的最优实验设计是一致的。这个性质在实际应用中非常有用，因为我们可以自由选择最方便的参数表示形式而不影响最优设计的结果。
在D-最优设计的实践中，最常见的情况是对少数几种不同的实验条件进行重复测量，正如案例3所展示的那样。这种重复测量的策略能够在有限的实验条件选择下，通过增加关键条件下的测量次数来最大化信息增益。
Fisher信息矩阵的构造
基本推导过程
考虑观测模型 y(i) = y_m(\mathbf{p}, i) + \varepsilon_i，其中 y_m(\mathbf{p}, i) 是由参数 \mathbf{p} 和实验条件 i 决定的确定性模型预测值，\varepsilon_i 是测量误差，服从独立同分布的概率密度 \pi(\varepsilon)。根据Fisher信息矩阵的定义，我们需要计算对数似然函数的二阶偏导：
由于观测值 y(i) 可以写成 y(i) = y_m(\mathbf{p}, i) + \varepsilon_i，误差项的概率密度为 \pi(\varepsilon_i)，因此观测值的概率密度为 \pi[y(i) - y_m(\mathbf{p}, i)]。对数似然函数对参数的导数可以通过链式法则展开。经过推导，可以将期望中的积分表达式展开为：
这个表达式中，积分项 \int \left(\frac{d\log \pi(e)}{de}\right)^2 \pi(e)de 是关于误差分布 \pi 的一个常数，它刻画了误差分布本身包含的信息量，这个量称为误差分布的Fisher信息。模型对参数的偏导数 \frac{\partial y_m(\mathbf{p}, i)}{\partial \mathbf{p}} 描述了当参数变化时模型预测值如何变化，它反映了模型对参数的敏感程度。
因此，Fisher信息矩阵可以简洁地表示为：
这里引入了两个关键概念。敏感度函数 \mathbf{s}_y(\mathbf{p}, i) 定义为模型输出对参数的梯度 \frac{\partial y_m(\mathbf{p}, i)}{\partial \mathbf{p}}，它是一个列向量，每个分量表示模型输出对相应参数的偏导数。敏感度函数衡量了在实验条件 i 下，模型对各个参数的敏感程度。分布的Fisher信息 \mathcal{I}(\pi) 定义为：
这个量完全由误差分布决定，与模型和参数无关。对于正态分布 \pi = \mathcal{N}(0, \sigma^2)，可以直接计算得到 \mathcal{I}(\pi) = 1/\sigma^2。这个结果说明，测量误差的方差越小，每次观测包含的信息量越大。
非平稳噪声的情形
如果测量误差 \varepsilon_i 不是平稳的，即不同观测的误差分布可能不同，那么每个观测 i 对应的误差有自己的分布 \pi_i。在这种情况下，Fisher信息矩阵的表达式变为：
这里 \pi_i 是第 i 次观测对应的误差概率密度。这个公式表明，不同观测对信息矩阵的贡献由两个因素决定：该观测条件下的敏感度 \mathbf{s}_y(\mathbf{p}, i) 和该观测的误差分布的Fisher信息 \mathcal{I}(\pi_i)。
但在后续的讨论中，我们假设误差是平稳的且服从正态分布，这个假设不失一般性。在平稳正态噪声假设下，所有观测的误差分布相同，因此 \mathcal{I}(\pi_i) 对所有 i 都等于 1/\sigma^2。这个假设简化了分析，使得Fisher信息矩阵只依赖于敏感度函数：
在这种情况下，Fisher信息矩阵正比于 1/\sigma^2，而优化问题中的目标函数（如行列式）只依赖于信息矩阵的相对大小，因此 \sigma^2 的具体数值在优化过程中可以约去，不影响最优设计的选择。
案例回顾与扩展
案例1的数学形式化
回到称量三个物体的问题，现在用统一的数学框架来表述。观测模型可以写成线性形式：
其中参数向量为 \mathbf{p}^* = (m_0, m_1, m_2, m_3)^{\mathrm{T}}，包含系统偏差 m_0 和三个物体的真实质量。设计矩阵 \mathbf{R} 描述了每次称量中各个参数的系数：
矩阵的第一列全为1，对应系统偏差 m_0 在每次称量中都存在。后三列对应三个物体的质量，r_{ij} = 1 表示第 i 次称量时物体 j 在天平上，r_{ij} = 0 表示物体 j 不参与第 i 次称量。由于每个物体要么在盘上要么不在，设计矩阵的元素只能取两个可能的值：0或1（对于单盘天平），或者-1、0、1（对于双盘天平）。
根据线性模型的Fisher信息矩阵公式，对于平稳正态噪声，信息矩阵为：
这里 N 是称量次数。对于单盘天平，系统偏差 m_0 总是存在的，无法通过单独的称量消除，因此可以将 m_0 视为干扰参数。如果我们只关心三个物体质量的估计精度，可以使用 D_s-最优性准则。
第一个注记指出，可以将 m_0 作为干扰参数处理，使用 D_s-最优性来设计实验。第二个注记强调，设计矩阵 \mathbf{R} 的元素只有两个可能取值（0或1），这使得这个问题具有非常特殊的组合优化结构。这是一个典型的离散优化问题。
在后续讨论中，实验变量 \boldsymbol{\xi} 将采用更一般的定量取值形式，不再局限于0和1这样的离散值，而是可以在连续区间内取值。
案例4：指数衰减模型
考虑一个简单的指数衰减过程，模型为：
其中 \varepsilon_i 是独立同分布的正态误差，\varepsilon_i \sim \mathcal{N}(0, \sigma^2)。这个模型有两个参数：\mathbf{p}_1 是初始幅度，\mathbf{p}_2 是衰减率。假设只进行两次观测，实验设计变量为 \boldsymbol{\xi}_1^2 = (t_1, t_2)，其中 0 \leq t_1 \leq t_2，即第一次观测时间不晚于第二次。
模型对参数的敏感度函数为：
第一个分量 \exp(-\mathbf{p}_2 t) 是模型对 \mathbf{p}_1 的偏导数，第二个分量 -t\mathbf{p}_1 \exp(-\mathbf{p}_2 t) 是模型对 \mathbf{p}_2 的偏导数。将两次观测的敏感度函数代入Fisher信息矩阵公式，经过计算可以得到：
行列式的表达式包含三项的乘积：\mathbf{p}_1^2 反映了信号幅度的影响，(t_2 - t_1)^2 说明两次观测时间的间隔越大信息量越大，指数项 \exp[-2\mathbf{p}_2(t_1 + t_2)] 表明观测时间的总和影响信息量。
为了最大化这个行列式，需要选择最优的观测时间。对行列式关于 t_1 和 t_2 求导并令其为零，可以解得D-最优设计为：
最优设计是在时间0（过程开始时）进行第一次观测，在 t = 1/\mathbf{p}_2（衰减时间常数的倒数）进行第二次观测。这个结果符合直觉：初始时刻信号最强，能提供关于幅度 \mathbf{p}_1 的最多信息；在一个时间常数处观测能最好地估计衰减率 \mathbf{p}_2。
这个案例揭示了实验设计的一个根本困难：最优设计 \boldsymbol{\xi}_D 依赖于真实参数 \mathbf{p}^* 的值，而参数恰恰是我们需要通过实验来估计的未知量。这种循环依赖在非线性模型的实验设计中总是存在。
精确实验设计
问题的基本设定
精确实验设计处理的是以下优化问题：对于 N 次观测 y(i)，需要为每次观测选择实验条件 \xi_i（称为第 i 个支撑点）。实验设计的目标是从初始设计 (\xi_1, \ldots, \xi_N) 优化到最优设计 \boldsymbol{\xi}_1^N，使得每个 \xi_i 都被单独优化。
这个问题的维度定义为 d = N \times \dim(\xi_i)，即问题的总变量数等于观测次数乘以每个实验条件的维度。例如，如果每个实验条件是一个时间点（一维），进行10次观测，那么 d = 10；如果每个实验条件包含温度和压力两个变量（二维），进行10次观测，那么 d = 20。
优化算法的选择取决于问题的维度。当 d 不太大时（比如 d < 100），可以使用经典的优化算法，即使存在约束条件和局部极小值问题，这些算法仍然可以有效处理。约束条件通常来自实验条件的物理限制，例如时间必须非负、温度有上下界等。局部极小值问题在非凸优化中普遍存在，需要使用多起点优化或全局优化方法来缓解。
当 d 非常大时，经典优化算法会遇到计算困难，需要开发专门的算法。这些专门算法通常利用问题的特殊结构，比如信息矩阵的可加性、敏感度函数的特殊形式等，来提高计算效率。大规模实验设计问题在工业应用中很常见，特别是在需要大量观测的复杂系统中。
交换算法的基本思想
交换算法是求解精确实验设计问题的一类迭代优化方法。算法的核心思想是在每次迭代中，从当前设计 \boldsymbol{\xi}^k = (\xi_1, \ldots, \xi_i, \ldots, \xi_N) 出发，选择其中一个支撑点 \xi_i 进行替换。替换的原则是找到一个更好的实验条件 \xi_j^*，使得替换后的新设计在准则函数 \Phi(\mathbf{M}) 的意义下优于原设计。这里"更好"是指准则函数值更小（对于需要最小化的准则）或更大（对于需要最大化的准则，如D-最优性中的行列式）。
算法示意如下：在第 k 次迭代时，当前设计为 \boldsymbol{\xi}^k = (\xi_1, \ldots, \xi_i, \ldots, \xi_N)，我们选择其中的 \xi_i 并用 \xi_j^* 替换它，得到新的设计。这个过程可以表示为：
箭头表示 \xi_i 被替换为 \xi_j^*。关键问题是如何选择要替换的支撑点 \xi_i 和新的实验条件 \xi_j^*。
Fedorov交换算法
Fedorov在1972年提出的算法采用穷举搜索策略。在每次迭代中，算法考虑所有可能的 N 种交换方式：依次尝试替换 \xi_1, \xi_2, \ldots, \xi_N 中的每一个支撑点。对于每个支撑点 \xi_j，算法寻找最优的替换实验条件 \xi_j^*，使得准则函数最优。这个优化过程可以表示为：
每个箭头对应一次维度为 \dim(\xi_i) 的优化问题。算法对每个位置 j = 1, \ldots, N 都进行这样的优化，得到 N 个候选替换方案。然后比较这 N 个替换方案对应的准则函数值，选择其中改进最大的那个作为本次迭代的实际替换，从而得到 \boldsymbol{\xi}^{k+1}。
Fedorov算法的特点是每次迭代需要进行 N 次维度为 \dim(\xi_i) 的优化，然后在这 N 个优化结果中比较 N 个准则函数值。这种穷举式的搜索保证了算法在当前设计的邻域内找到最优的单点替换，但计算代价较高。当观测次数 N 较大时，每次迭代的计算负担会变得很重。
DETMAX算法
Mitchell在1974年提出的DETMAX算法采用了一种更高效的策略。算法的基本思路是通过增加和删除支撑点的组合来寻找改进方向。具体来说，算法假设可以进行一次额外的实验，即暂时将设计从 N 个观测扩展到 N+1 个观测。算法首先对第 N+1 个支撑点 \xi_{N+1} 进行优化，找到最优的 \xi_{N+1}^*：
这里 \xi_{N+1} 通过优化确定，使得包含 N+1 个观测的设计达到最优。这只需要一次维度为 \dim(\xi_i) 的优化。
接下来，算法需要回到 N 个观测的约束，因此必须删除一个支撑点。算法考虑所有 N+1 种可能的删除方案，即依次尝试删除 \xi_1, \xi_2, \ldots, \xi_N, \xi_{N+1} 中的每一个。对于每种删除方案，计算删除后的准则函数值。算法选择删除后对准则函数影响最小的那个支撑点，设为 \xi_j。这相当于全局地用新加入的 \xi_{N+1}^* 交换掉原设计中的 \xi_j。
这种先增后减的策略构成了长度为1的游走过程。在实践中，可以考虑更长的游走过程，即连续进行多次增加和删除操作，但这会增加计算复杂度。DETMAX算法的一次迭代只需要1次维度为 \dim(\xi_i) 的优化，加上 N+1 次准则函数值的比较。相比Fedorov算法，DETMAX的单次迭代计算量更小，因此迭代速度更快。但是，DETMAX通常需要更多的迭代次数才能收敛到最优设计。
算法的死锁问题
交换算法可能会遇到死锁现象，即算法陷入局部最优而无法继续改进。对于Fedorov算法，死锁发生的条件是无法通过优化任何单个支撑点 \xi_i 来改进准则函数。换句话说，对于当前设计中的每个 \xi_i，任何替换都不能使准则函数变得更好。
对于DETMAX算法，死锁的条件是要删除的支撑点恰好是刚刚添加的 \xi_{N+1}。这意味着算法找到的最优新支撑点 \xi_{N+1}^* 对设计没有贡献，删除它对准则函数的影响最小。此时算法退化为原地踏步，无法产生新的设计。
死锁现象反映了优化问题的非凸性和局部极小值的存在。当算法收敛到局部最优时，需要使用多起点策略或其他全局优化技术来跳出局部最优，寻找更好的设计。
近似实验设计
基本概念和定义
近似实验设计是精确实验设计的一种松弛形式。在精确设计中，Fisher信息矩阵的完整表达式为：
这个求和包含 N 项，每项对应一次观测。当多个观测使用相同的实验条件时，即存在重复测量，可以将具有相同实验条件的观测项合并。假设在 m 个不同的实验条件 \xi_i 下进行测量，其中 m < N，每个条件 \xi_i 下进行 r_i 次重复测量，满足 \sum_{i=1}^{m} r_i = N。此时信息矩阵可以重新表示为：
这里 \frac{r_i}{N} 表示在实验条件 \xi_i 下收集的观测比例。定义 \lambda(\xi_i) = \frac{r_i}{N} 为实验条件 \xi_i 的权重，它表示在该条件下进行测量的相对频率。权重满足归一化条件 \sum_{i=1}^{m} \lambda(\xi_i) = 1 和非负性约束 \lambda(\xi_i) \geq 0。这个权重函数 \lambda(\xi_i) 在后续讨论中将被记作离散概率测度。
近似实验设计将权重从离散的有理数 \frac{r_i}{N} 推广到连续的实数 \lambda(\xi_i) \in [0, 1]。这种推广允许权重取任意实数值，不再受整数观测次数的限制。近似设计的优势在于它将一个困难的整数规划问题（精确设计）转化为一个相对容易的连续优化问题，可以利用凸优化等强大的数学工具来求解。
近似设计的表示形式
在近似实验设计框架下，平均每次观测的Fisher信息矩阵表示为：
实验方案 \boldsymbol{\xi}_1^N 现在可以用一个紧凑的形式表示：
上面一行列出所有不同的实验条件（支撑点），下面一行给出每个条件对应的权重。权重满足归一化条件 \sum_{i=1}^{m} \lambda(\xi_i) = 1，这形成了实验条件 \xi_i 上的一个离散归一化分布。由于在精确设计中 \lambda(\xi_i) = \frac{r_i}{N} 必须是有理数，这给优化带来了整数约束的困难。
近似设计的核心思想是进行连续松弛，即允许权重 \lambda_i 取任意非负实数，只要满足：
其中 0 \leq \lambda_i，i = 1, \ldots, m，且 \sum_{i=1}^{m} \lambda_i = 1。这样的 \boldsymbol{\xi} 被称为在实验条件空间上的离散归一化分布，也可以理解为一个离散概率分布，其中支撑点为 \xi_i，权重为 \lambda_i。在概率论的语言中，这被称为统计方案。
从离散到连续设计
约束条件可以进一步放松，将 \boldsymbol{\xi} 从离散分布推广到连续的归一化测度。设 \mathcal{X} 为所有可行实验条件构成的集合（可行域），\boldsymbol{\xi} 是定义在 \mathcal{X} 上的归一化测度，不一定是离散的。此时Fisher信息矩阵表示为积分形式：
这里 \boldsymbol{\xi}(d\xi) 表示测度 \boldsymbol{\xi} 在微元 d\xi 上的测度值。积分条件为 \int_{\mathcal{X}} \boldsymbol{\xi}(d\xi) = 1，确保 \boldsymbol{\xi} 是归一化的。这种形式化来自Kiefer在1959年提出的连续协议概念，它将实验设计问题转化为在测度空间上的优化问题。
连续设计的引入极大地简化了数学分析，因为测度空间具有良好的凸性和拓扑性质，可以利用泛函分析和凸优化理论的强大工具。尽管实际实验中只能进行有限次离散的观测，连续设计提供了理论上的最优基准，离散设计可以看作是对连续最优设计的近似实现。
信息矩阵的凸包络结构
Fisher信息矩阵 \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 属于集合 \mathcal{M}_1 的凸包络。这里 \mathcal{M}_1 是所有秩为1的矩阵构成的集合，这些矩阵具有如下形式：
每个这样的矩阵对应于在单一实验条件 \boldsymbol{\xi} 下进行观测所获得的信息矩阵。由于敏感度函数 \mathbf{s}_y 是一个列向量，其外积 \mathbf{s}_y\mathbf{s}_y^{\mathrm{T}} 是一个秩为1的对称半正定矩阵。
任何通过多次观测获得的Fisher信息矩阵都可以表示为这些秩1矩阵的凸组合，即：
其中 \lambda_i \geq 0 且 \sum_{i=1}^{m} \lambda_i = 1。由于 \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 是对称矩阵，它属于维度为 d = \frac{n_p(n_p+1)}{2} 的向量空间，这里 n_p = \dim(\mathbf{p}) 是参数的个数。对称矩阵的独立元素个数等于上三角（包括对角线）的元素个数，即 \frac{n_p(n_p+1)}{2}。
图示展示了这个几何结构：\mathcal{M}_1 是所有秩1信息矩阵构成的集合，它通常是一个非凸的曲面或流形。凸包络 \text{env.conv}(\mathcal{M}_1) 是包含 \mathcal{M}_1 的最小凸集，所有可能的Fisher信息矩阵都落在这个凸包络内部或边界上。这个凸结构是近似实验设计理论的基础，它保证了优化问题具有良好的数学性质。
Carathéodory定理的应用
Carathéodory定理是凸分析中的一个基本结果。在 n 维仿射空间中，如果一个点属于集合 \mathcal{A} 的凸包络，那么这个点可以表示为 \mathcal{A} 中至多 n+1 个点的重心，即这些点的凸组合，且组合系数为非负或零。
将这个定理应用到实验设计问题，\mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 属于 \mathcal{M}_1 的凸包络，而 \mathcal{M}_1 位于维度为 d = \frac{n_p(n_p+1)}{2} 的空间中。因此，任何Fisher信息矩阵可以表示为至多 d+1 个秩1矩阵的线性组合：
这个结果意味着，在近似设计框架下，测度 \boldsymbol{\xi} 对应的是一个至多有 \frac{n_p(n_p+1)}{2} + 1 个支撑点的离散测度。即使我们允许测度是连续的，最优设计总是可以用一个有限个支撑点的离散测度来实现。
Carathéodory定理的一个重要推论是：这个结论对最优方案同样成立。无论我们使用什么样的凸优化准则（如D-最优性、A-最优性等），最优设计总是可以用至多 \frac{n_p(n_p+1)}{2} + 1 个不同的实验条件来实现。这大大简化了最优设计的搜索空间，从无限维的测度空间降低到有限维的组合优化问题。例如，对于两个参数的问题（n_p = 2），d = \frac{2 \times 3}{2} = 3，最优设计至多需要4个不同的实验条件。对于三个参数（n_p = 3），d = 6，最优设计至多需要7个支撑点。
D-最优设计的边界性质
现在进一步探讨D-最优方案的一个重要几何性质。设 \xi_D 是D-最优协议，即它最大化 \det \mathbf{M}(\mathbf{p}, \boldsymbol{\xi})。假设对应的Fisher信息矩阵 \mathbf{M}(\mathbf{p}, \xi_D) 位于 \mathcal{M}_1 的凸包络的内部。这意味着 \mathbf{M}(\mathbf{p}, \xi_D) 不在边界上，可以在其周围的任意方向上进行小的扰动而仍然保持在凸包络内。
如果 \mathbf{M}(\mathbf{p}, \xi_D) 位于内部，那么必然存在某个标量 a > 1，使得缩放后的矩阵 a\mathbf{M}(\mathbf{p}, \xi_D) 恰好落在 \mathcal{M}_1 的边界上。这是因为从内部任一点出发，沿着任意方向延伸，最终必然会到达边界。计算缩放矩阵的行列式：
这里利用了行列式的齐次性：对矩阵乘以标量 a，其行列式变为 a^{n_p}，其中 n_p 是矩阵的维度（参数个数）。由于 a > 1，缩放后的矩阵行列式严格大于原矩阵。但这产生了矛盾：如果存在另一个可达的信息矩阵 a\mathbf{M}(\mathbf{p}, \xi_D) 其行列式更大，那么 \xi_D 就不应该是D-最优方案。
这个矛盾说明初始假设不成立，因此D-最优设计对应的Fisher信息矩阵必然位于凸包络的边界上：
这个几何性质表明，最优设计总是"紧绷"在可行域的边界上，充分利用了可行信息矩阵集合的极限。
从近似设计到精确设计的实现
Carathéodory定理的延续告诉我们，D-最优信息矩阵 \mathbf{M}(\mathbf{p}, \xi_D) 可以表示为至多 d 个 \mathcal{M}_1 中元素的线性组合：
这给出了最优近似设计的支撑点数量上界。但近似设计中的权重 \lambda_i 可以是任意实数，而精确实验只能进行整数次观测。因此需要将近似设计转化为可实施的精确设计。
转化方法是考虑一个包含重复测量的精确协议。设总观测次数为 N，在实验条件 \xi_i 下进行 r_i 次重复测量，满足 \sum_{i=1}^{m} r_i = N。精确设计的权重为：
为了使精确设计接近最优近似设计，我们选择一个足够大的 N，然后选择整数 r_i 使得 \frac{r_i}{N} 尽可能接近最优权重 \lambda_i。当 N 趋于无穷时，可以让 \frac{r_i}{N} 任意接近 \lambda_i，这样得到的精确设计称为渐近方案。
这里需要注意，在某些特殊情况下，可以直接选择 \boldsymbol{\xi} 而无需进行任何近似。例如，当实验条件 \boldsymbol{\xi} 代表信号的平均功率谱密度时，测度 \boldsymbol{\xi} 本身就是一个连续的谱分布，可以直接在连续域上进行优化和实施。
近似实验设计的凸性质
近似实验设计框架具有优美的数学性质。首先，可行信息矩阵的集合 \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 构成一个凸集。这是因为对于任意两个设计 \boldsymbol{\xi}_1 和 \boldsymbol{\xi}_2，它们的凸组合：
仍然是一个有效的设计（一个归一化测度），对应的信息矩阵满足线性关系：
更重要的是，D-最优性的目标函数，即对数行列式 \log \det \mathbf{M}，是凹函数。对于任意两个正定矩阵 \mathbf{M}_1 和 \mathbf{M}_2，且 \mathbf{M}_1 \neq \mathbf{M}_2，对于任意 0 < \alpha < 1，可以证明：
不等式右边是两个对数行列式值的线性插值，左边是凸组合矩阵的对数行列式。凹函数的性质保证了线性插值总是低于函数值本身。这意味着 \log \det[\cdot] 是一个凹函数。
凸集上的凹函数优化问题具有非常好的性质：局部最优解就是全局最优解，不存在其他局部极大值陷阱。这为D-最优实验设计的数值求解提供了理论保证，使用梯度法、牛顿法等标准凸优化算法都能有效收敛到全局最优解。
D-最优性的充要条件推导
现在证明 \xi_D 是D-最优方案当且仅当它最大化 \log \det \mathbf{M}(\mathbf{p}, \boldsymbol{\xi})。考虑从 \xi_D 出发沿任意方向 \tilde{\boldsymbol{\xi}} 的扰动，构造凸组合：
其中 0 \leq \alpha \leq 1。由于信息矩阵关于设计的线性性，扰动后的信息矩阵为：
\xi_D 是D-最优的充要条件是在所有可能的扰动方向上，目标函数都不能进一步增加。用一阶最优性条件表述，即在 \alpha = 0 处，目标函数 \log \det \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 关于 \alpha 的导数必须非正：
这个条件的几何意义是：如果 \xi_D 确实是最优的，那么在其位置处，目标函数在所有方向上的斜率都不应该是正的，否则我们可以沿着正斜率方向移动来进一步改进。
D-最优性的几何理解
图示展示了这个优化问题的几何结构。横轴代表设计空间，纵轴表示 \log \det \mathbf{M}(\mathbf{p}, \boldsymbol{\xi}) 的值。可行信息矩阵的集合 \mathcal{M}_1 及其凸包络形成了一个曲面，D-最优设计 \mathbf{M}(\mathbf{p}, \xi_D) 对应于这个曲面上使得对数行列式达到最大值的点。
这里有一个重要的注记：D-最优信息矩阵 \mathbf{M}(\mathbf{p}, \xi_D) 是唯一的，这是因为目标函数 \log \det \mathbf{M} 是严格凹函数，严格凹函数在凸集上的最大值点是唯一的。但是，D-最优设计 \xi_D 本身不一定唯一。这看似矛盾，但实际上反映了一个事实：可能存在多个不同的测度 \xi_D 产生相同的信息矩阵。D-最优设计的集合本身构成一个凸集，即如果 \xi_D^{(1)} 和 \xi_D^{(2)} 都是D-最优的，那么它们的任意凸组合也是D-最优的。
"我们在顶部"的含义是：D-最优点对应于目标函数的全局最大值，在该点处函数达到峰值，任何方向的移动都会导致函数值下降或至多保持不变。
方向导数的计算
利用矩阵微分的链式法则，可以计算对数行列式关于参数 \alpha 的导数。对于可逆矩阵 \mathbf{M}，有恒等式：
这个公式来源于 \frac{d \log \det \mathbf{M}}{d\mathbf{M}} = \mathbf{M}^{-\mathrm{T}}。在 \alpha = 0 处，\boldsymbol{\xi} = \xi_D，而 \frac{\partial \mathbf{M}}{\partial \alpha}\Big|_{\alpha=0} = \mathbf{M}(\mathbf{p}, \tilde{\boldsymbol{\xi}}) - \mathbf{M}(\mathbf{p}, \xi_D)。代入得到：
利用迹的线性性和 \text{trace}(\mathbf{M}^{-1}\mathbf{M}) = n_p（因为 \mathbf{M}^{-1}\mathbf{M} = \mathbf{I}，单位矩阵的迹等于其维度），条件可以简化为：
这个条件必须对所有可能的设计 \tilde{\boldsymbol{\xi}} 都成立。特别地，考虑退化情况：选择 \tilde{\boldsymbol{\xi}} = \delta_{\xi}，即在单个实验条件 \xi 处集中全部测量的狄拉克测度。此时：
这是一个秩为1的矩阵。代入得到D-最优性的必要条件：
这里利用了循环迹的性质 \text{trace}(\mathbf{ABC}) = \text{trace}(\mathbf{BCA}) 和 \text{trace}(\mathbf{vv}^{\mathrm{T}}) = \mathbf{v}^{\mathrm{T}}\mathbf{v}。
充分性证明
现在证明上述条件也是充分的。假设对所有 \xi \in \mathcal{X}，都有：
对于任意设计 \tilde{\boldsymbol{\xi}}，根据Carathéodory定理，它可以表示为至多 \tilde{m} \leq \frac{n_p(n_p+1)}{2} + 1 个支撑点的凸组合：
计算迹：
最后一步利用了 \sum_{i=1}^{\tilde{m}} \lambda_i = 1。因此，如果点条件对所有 \xi \in \mathcal{X} 成立，那么对所有设计 \tilde{\boldsymbol{\xi}} 也成立，这证明了 \xi_D 是D-最优的。
总结：D-最优性的充要条件为，对所有可行实验条件 \xi \in \mathcal{X}，定义灵敏度函数：
如果某个设计 \xi 不是D-最优的，那么必然存在某个 \xi^* \in \mathcal{X} 使得 d(\xi^*, \xi) > n_p。
Kiefer-Wolfowitz等价定理
Kiefer和Wolfowitz在1960年证明了一个优美的等价定理，它将D-最优性用三种不同但等价的方式刻画：
(1) \xi_D 是D-最优的
(2) \max_{\xi \in \mathcal{X}} d(\xi, \xi_D) = n_p
(3) \xi_D 最小化 \max_{\xi \in \mathcal{X}} d(\xi, \xi_D)
这三个条件相互等价：(1) \Leftrightarrow (2) \Leftrightarrow (3)。
第二个条件说明，在D-最优设计处，灵敏度函数 d(\xi, \xi_D) 在整个可行域上的最大值恰好等于参数个数 n_p。这个最大值在 \xi_D 的所有支撑点处达到。事实上，对于 \xi_D 的任意支撑点 \xi_i（即权重 \lambda_i > 0 的点），都有 d(\xi_i, \xi_D) = n_p。这可以通过计算验证：
由于 d(\xi_i, \xi_D) \leq n_p 对所有 i 成立，且加权平均等于 n_p，所以每个支撑点处必有 d(\xi_i, \xi_D) = n_p。
第三个条件给出了D-最优设计的极小极大解释：它在所有可能的设计中，最小化了灵敏度函数的最大值。这个解释与鲁棒性相关，说明D-最优设计使得在最坏情况下的预测不确定性达到最小。
统计意义
等价定理还揭示了D-最优性与预测方差的关系。如果用协议 \boldsymbol{\xi}_1^N 估计参数得到 \hat{\mathbf{p}}，那么在新实验条件 \xi 处的模型预测 y_m(\hat{\mathbf{p}}, \xi) 的方差可以近似为：
利用泰勒展开和Fisher信息矩阵作为渐近协方差矩阵的逆，可以得到：
注意到 \frac{\partial y_m(\mathbf{p}, \xi)}{\partial \mathbf{p}} = \mathbf{s}_y(\mathbf{p}, \xi)，因此预测方差正比于 d(\xi, \boldsymbol{\xi}_1^N)：
这意味着条件(3)可以重新解释为：\xi_D 最小化了在所有可能实验条件 \xi \in \mathcal{X} 上，模型输出预测方差的最大值。D-最优设计保证了在最不利的预测位置，不确定性也被控制在最小水平，这体现了设计的鲁棒性和均衡性。