企业官网建设流程全解析

• 科学计算

导读

本文基于 CVXPY 官方示例 queuing_design.rst，完整讲解如何将 M/M/N 排队系统的服务负载最小化问题建模为 log-log 凸规划（LLCP），并使用 CVXPY 的requires_grad=True可微求解能力，通过problem.backward()与problem.derivative()分别做参数灵敏度分析与一阶扰动预测。读完本文，你将掌握 DGP/LLCP 建模、backward/derivative两种微分模式的 API 语义，以及如何用源码级机制（DPP 合规检查、diffcp 求解器缓存、链式法则约简）验证分析结果的可靠性。

M/M/N 排队系统：问题背景与数学模型

排队系统是一组队列的集合，队列中的任务等待被服务——这些任务可能是操作系统中的线程，也可能是网络系统输入/输出缓冲区中的数据包。本文考虑一个包含 $N$ 个队列的马尔可夫排队系统（M/M/N 队列），设计目标是在延迟、占用等约束下最小化系统的服务负载。

假设到达第 $i$ 个队列的任务由速率为 $\lambda_i$ 的 Poisson 过程产生，第 $i$ 个队列的服务时间服从参数为 $\mu_i$ 的指数分布，$i=1,\ldots,N$。系统核心度量是服务负载函数 $\ell: \mathbf{R}^{N}{++} \times \mathbf{R}^{N}{++} \to \mathbf{R}^{N}_{++}$：

$$ \ell_i(\lambda, \mu) = \frac{\mu_i}{\lambda_i}, \quad i=1, \ldots, N. $$

注意这是交通负荷 $\rho$（通常记作 $\rho$）的倒数。在此基础上，系统还定义了三组关键性能函数：队列占用率 $q$、平均延迟 $w$ 与总延迟 $d$：

$$ q_i(\lambda, \mu) = \frac{\ell_i(\lambda, \mu)^{-2}}{1 - \ell_i(\lambda, \mu)^{-1}}, \quad w_i(\lambda, \mu) = \frac{q_i(\lambda, \mu)}{\lambda_i} + \frac{1}{\mu_i}, \quad d_i(\lambda, \mu) = \frac{1}{\mu_i - \lambda_i} $$

这些函数的定义域为 ${(\lambda, \mu) \in \mathbf{R}^{N}{++} \times \mathbf{R}^{N}{++} \mid \lambda < \mu }$（逐元素不等式），即到达率必须严格小于服务率，这是排队系统稳定（不无限堆积）的基本前提。

系统受到的约束为：队列占用率、平均排队延迟与总延迟分别不超过参数 $q_{\max}$、$w_{\max}$、$d_{\max}$（均为 $\mathbf{R}^{N}{++}$ 元素，不等式逐元素成立）；同时到达率向量 $\lambda$ 至少为 $\lambda{\mathrm{min}} \in \mathbf{R}^{N}{++}$，且服务率总和不超过 $\mu{\max} \in \mathbf{R}_{++}$。

将设计问题建模为 Log-Log 凸规划（LLCP）

设计问题为：在满足上述约束的前提下，最小化服务负载的加权和 $\gamma^T \ell(\lambda, \mu)$，其中 $\gamma \in \mathbf{R}^{N}_{++}$ 是权重向量：

$$ \begin{array}{ll} \mbox{minimize} & \gamma^T \ell(\lambda, \mu) \ \mbox{subject to} & q(\lambda, \mu) \leq q_{\max} \ & w(\lambda, \mu) \leq w_{\max} \ & d(\lambda, \mu) \leq d_{\max} \ & \lambda \geq \lambda_{\mathrm{min}}, \quad \sum_{i=1}^{N} \mu_i \leq \mu_{\max}. \end{array} $$

其中 $\lambda, \mu \in \mathbf{R}^{N}{++}$ 是优化变量；$\gamma, q{\max}, w_{\max}, d_{\max}, \lambda_{\mathrm{min}} \in \mathbf{R}^{N}{++}$ 以及 $\mu{\max} \in \mathbf{R}_{++}$ 是参数。

该问题是一个 LLCP。目标函数 $\gamma^T \ell(\lambda, \mu)$ 是 posynomial（正项式），约束函数 $w$ 同样是 posynomial。而 $d$ 和 $q$ 不是 posynomial，但它们是log-log 凸的：

• $d$ 的 log-log 凸性由复合规则得出：函数 $(x, y) \mapsto y - x$ 在 $0 < x < y$ 上 log-log 凹，比值 $(x, y) \mapsto x/y$ 是 log-log 仿射且关于 $y$ 递减，因此 $d_i = 1/(\mu_i - \lambda_i)$ 是 log-log 凸的；
• 类似地，$q$ 也是 log-log 凸的。

这一数学性质在 CVXPY 中由具体的原子实现保证。查看源码 cvxpy/atoms/one_minus_pos.py，diff_pos(x, y)被定义为multiply(x, one_minus_pos(y/x))，其 docstring 明确声明"该原子是 log-log 凹的"；而 one_minus_pos 类的is_atom_log_log_convex()返回False、is_atom_log_log_concave()返回True。模型中的d = 1/cp.diff_pos(mu, lam)与lq = (ell)**(-2)/cp.one_minus_pos(ell**(-1))正是建立在这些原子的 log-log 凹凸性之上的——DGP 的符号规则（log-log 凸性的复合规则）会据此判定整个表达式的 log-log 凸性。

CVXPY 实现：变量、参数与约束的完整代码

以下代码完整继承自官方示例，可直接复制运行。首先导入依赖并定义变量与参数：

import cvxpy as cp import numpy as np import time mu = cp.Variable(pos=True, shape=(2,), name='mu') lam = cp.Variable(pos=True, shape=(2,), name='lambda') ell = cp.Variable(pos=True, shape=(2,), name='ell') w_max = cp.Parameter(pos=True, shape=(2,), value=np.array([2.5, 3.0]), name='w_max') d_max = cp.Parameter(pos=True, shape=(2,), value=np.array([2., 2.]), name='d_max') q_max = cp.Parameter(pos=True, shape=(2,), value=np.array([4., 5.0]), name='q_max') lam_min = cp.Parameter(pos=True, shape=(2,), value=np.array([0.5, 0.8]), name='lambda_min') mu_max = cp.Parameter(pos=True, value=3.0, name='mu_max') gamma = cp.Parameter(pos=True, shape=(2,), value=np.array([1.0, 2.0]), name='gamma')

这里的关键点：

• 所有变量均以pos=True声明，强制取正，对应数学定义域 $\mathbf{R}^{N}_{++}$；
• 所有参数也以pos=True声明，同样约束为正；向量参数使用shape=(2,)（本示例取 $N=2$ 个队列），标量参数mu_max与gamma分别控制服务率总预算与目标加权；
• 给参数显式name是为了后续打印梯度时能清晰识别每个参数的输出。

接着用原子组合出性能函数与约束：

lq = (ell)**(-2)/cp.one_minus_pos(ell**(-1)) q = lq w = lq/lam + 1/mu d = 1/cp.diff_pos(mu, lam) constraints = [ w <= w_max, d <= d_max, q <= q_max, lam >= lam_min, cp.sum(mu) <= mu_max, ell == mu/lam, ] objective_fn = gamma.T @ ell problem = cp.Problem(cp.Minimize(objective_fn), constraints)

对应关系一目了然：lq即 $q_i = \ell_i^{-2}/(1-\ell_i^{-1})$，w = lq/lam + 1/mu即 $w_i = q_i/\lambda_i + 1/\mu_i$，d = 1/cp.diff_pos(mu, lam)即 $d_i = 1/(\mu_i-\lambda_i)$。约束列表完整对应数学建模中的五组不等式，其中ell == mu/lam是服务负载的定义等式，把 $\ell$ 作为显式变量引入，便于后续对其单独做灵敏度分析。

求解：以 DGP 模式开启梯度支持

problem.solve(requires_grad=True, gp=True, eps=1e-14, max_iters=10000, mode='dense')

求解结果（目标函数最优值）为：

4.457106781186705

最优解打印如下：

print('mu ', mu.value) print('lam ', lam.value) print('ell ', ell.value)

mu [1.32842713 1.67157287] lam [0.82842712 1.17157287] ell [1.60355339 1.4267767 ]

注意solve(requires_grad=True, gp=True, ...)的求解器调用细节：gp=True表示按 DGP/LLCP 模式解析问题，requires_grad=True则要求在求解时缓存可微求解器的内部数据，供后续backward()/derivative()使用。从 cvxpy/problems/problem.py 的源码可以看到启用requires_grad的三条硬性约束：

1. 问题必须是DPP 合规的（problem.is_dpp(dpp_context)，dpp_context在gp=True时为'dgp'），否则抛出DPPError；
2. 求解器必须为SCS 或 DIFFCP，否则抛出ValueError；
3. 若选用DIFFCP，需要安装独立的diffcpPython 包。

eps=1e-14与max_iters=10000是透传给 SCS 的收敛容差与迭代上限，mode='dense'指定使用稠密后端。从实现上看，requires_grad=True的求解结果会被缓存到self._solver_cache[s.DIFFCP]中（backward_cache），其中保存了求解器导出的微分算子D/DT，这正是后续两类微分 API 的数据基础。

灵敏度分析：用backward()计算变量对参数的梯度

在 LLCP 设计问题中，我们关心的核心问题是：最优设计对哪些参数最敏感？例如，放宽最大总延迟 $d_{\max}$ 或提高服务率总预算 $\mu_{\max}$ 对服务负载分别有多大影响？这可以通过backward()求解。

problem.solve(requires_grad=True, gp=True, eps=1e-14, max_iters=10000, mode='dense') for var in [lam, mu, ell]: print('Variable ', var.name()) print('Gradient with respect to first component') var.gradient = np.array([1., 0.]) problem.backward() for param in problem.parameters(): if np.prod(param.shape) == 2: print('{0}: {1:.3g}, {2:.3g}'.format(param.name(), param.gradient[0], param.gradient[1])) else: print('{0}: {1:.3g}'.format(param.name(), param.gradient)) print('Gradient with respect to second component') var.gradient = np.array([0., 1.]) problem.backward() for param in problem.parameters(): if np.prod(param.shape) == 2: print('{0}: {1:.3g}, {2:.3g}'.format(param.name(), param.gradient[0], param.gradient[1])) else: print('{0}: {1:.3g}'.format(param.name(), param.gradient)) var.gradient = np.zeros(2) print('')

运行输出（完整继承原文结果）：

Variable lambda Gradient with respect to first component gamma: 0.213, -0.107 w_max: 5.43e-12, 5.64e-12 d_max: -0.411, -0.113 q_max: 5.99e-12, 4.77e-12 lambda_min: -1.56e-11, -7.35e-12 mu_max: 0.927 Gradient with respect to second component gamma: -0.458, 0.229 w_max: 2.08e-11, 2.16e-11 d_max: -0.105, -0.466 q_max: 2.29e-11, 1.83e-11 lambda_min: -5.97e-11, -2.82e-11 mu_max: 1.01 Variable mu Gradient with respect to first component gamma: 0.213, -0.107 w_max: 1.55e-11, 1.6e-11 d_max: -0.661, -0.113 q_max: 1.7e-11, 1.36e-11 lambda_min: -4.43e-11, -2.09e-11 mu_max: -0.0727 Gradient with respect to second component gamma: -0.458, 0.229 w_max: 2.3e-11, 2.39e-11 d_max: -0.105, -0.716 q_max: 2.53e-11, 2.02e-11 lambda_min: -6.59e-11, -3.11e-11 mu_max: 0.00996 Variable ell Gradient with respect to first component gamma: -0.245, 0.122 w_max: 2e-11, 2.08e-11 d_max: -0.282, -0.22 q_max: 2.21e-11, 1.76e-11 lambda_min: -5.74e-11, -2.71e-11 mu_max: -0.334 Gradient with respect to second component gamma: 0.122, -0.0611 w_max: -1.24e-13, -1.29e-13 d_max: -0.101, -0.195 q_max: -1.37e-13, -1.09e-13 lambda_min: 3.58e-13, 1.66e-13 mu_max: -0.197

结果解读（原文的核心结论）：解对 $d_{\max}$ 和 $\mu_{\max}$ 高度敏感。以lambda为例，$\partial \lambda_1 / \partial d_{\max}$ 为-0.411、$\partial \lambda_1 / \partial \mu_{\max}$ 为0.927，说明增大 $d_{\max}$（放宽延迟上限）会降低最优到达率，而增大 $\mu_{\max}$（提高服务率总预算）会提升到达率——且对第一个队列的影响尤其显著。相比之下，$w_{\max}$、$q_{\max}$、$\lambda_{\min}$ 对应的梯度均为 $10^{-11}$ 量级，说明这些约束在当前最优解处非活跃（不 binding），参数微扰几乎不影响解。

源码级机制：problem.backward() 计算"解的梯度关于参数的梯度"，即对变量梯度向量 $dz/dx$ 通过链式法则反传得到 $dz/dp$。其实现要点为：

• 设置变量gradient属性即为指定 $dz/dx$；未设置时默认广播为全 1 向量（对应 $f$ 取求和函数）；
• 沿solving_chain.reductions依次调用各约简的var_backward，将外层变量梯度传导到参数化问题内部；随后利用缓存的DT算子与apply_param_jac得到参数梯度，再沿逆序的约简链param_backward还原到原始参数空间；
• 副作用是填充每个Parameter的gradient属性。注意它只能在solve(requires_grad=True, ...)之后调用，且要求问题状态为有解（否则抛出SolverError）。

扰动分析：用derivative()做一阶预测并验证

灵敏度分析回答的是"参数变化一点点，解怎么变"的线性化问题；扰动分析则进一步用一阶近似预测参数扰动后的新解，并与重新求解的真实解对比，验证线性近似的精度。

problem.solve(requires_grad=True, gp=True, eps=1e-14, max_iters=10000, mode='dense') mu_value = mu.value lam_value = lam.value delta = 0.01 for param in problem.parameters(): param.delta = param.value * delta problem.derivative() lam_pred = (lam.delta / lam_value) * 100 mu_pred = (mu.delta / mu_value) * 100 print('lam predicted (percent change): ', lam_pred) print('mu predicted (percent change): ', mu_pred) for param in problem.parameters(): param._old_value = param.value param.value += param.delta problem.solve(cp.SCS, gp=True, eps=1e-14, max_iters=10000) lam_actual = ((lam.value - lam_value) / lam_value) * 100 mu_actual = ((mu.value - mu_value) / mu_value) * 100 print('lam actual (percent change): ', lam_actual) print('mu actual (percent change): ', mu_actual)

运行输出：

lam predicted (percent change): [2.32203282 1.77228841] mu predicted (percent change): [1.07166961 0.94304296] lam actual (percent change): [1.99504983 1.99504965] mu actual (percent change): [0.87148458 1.10213353]

这里的操作协议是：

1. 给每个参数设置param.delta = param.value * delta（即相对扰动 $1%$）；
2. 调用problem.derivative()，它会填充变量delta属性——即一阶预测的变量变化量；
3. 把扰动真正施加到参数上（param.value += param.delta），重新求解原始问题得到真实解；
4. 分别计算预测与实际的百分比变化并对比。

对比可见：一阶预测与真实变化的量级完全一致（约 $1%\sim 2%$ 量级），方向相同、数值接近，验证了灵敏度分析结论的可靠性；当然在 $1%$ 的有限扰动下也存在可观测的非线性偏差，这正是线性近似的固有误差。

源码级机制：problem.derivative() 是backward()的"正向模式"对应物——它把参数扰动 $\delta p$ 通过解映射的导数映射为变量最优值的扰动 $\delta x$。实现要点：

• 未设置delta的参数默认扰动为 0（np.broadcast_to(0.0, param.shape)）；
• 沿约简链正向调用param_forward把参数扰动传导到参数化问题，利用缓存的D算子计算内部解扰动，再经逆序的var_forward还原为外层变量delta；
• 副作用是填充每个Variable的delta属性；同样要求先以requires_grad=True求解且问题有解。

这两类 API 的互补语义可以总结为：backward()= 给定变量的"梯度种子"，反传出每个参数的梯度（用于灵敏度分析 / 集成进自动微分工具链）；derivative()= 给定参数的"扰动种子"，前传预测每个变量的最优值变化（用于 what-if 分析 / 一阶近似）。从 solve 的 docstring 可以看出，二者都要求问题为 DPP 合规且求解器为 SCS/diffcp，同时明确"DQCP 问题不支持梯度计算"。

从源码理解：DGP 约简如何支撑可微 LLCP

本示例能同时做到"按 DGP 语义建模"与"求导"，得益于 CVXPY 的两层机制：

• DGP → DCP 约简：cvxpy/reductions/dgp2dcp/dgp2dcp.py 中的Dgp2Dcp类将 DGP 问题约简为等价的 DCP 问题——它通过对数变量替换（每个变量 $x$ 替换为 $u = \log x$）把 posynomial 约束变成凸约束，并在求解后通过invert把对数空间的解映射回原始变量。该约简同时记录_id_to_var等映射关系，作为后续梯度反传时还原变量身份的inverse_data；
• 可微求解器缓存：requires_grad=True时，diffcp（或 SCS 的可微接口）在求解过程中保存灵敏度算子 $D$（正向）与 $DT$（反向），backward()与derivative()正是利用这些算子配合约简链上的var_backward/param_backward/param_forward/var_forward实现端到端微分。这也是为什么启用梯度时求解器被强制限定为 SCS/DIFFCP、且问题必须满足 DPP。

更细节地，本示例用到的diff_pos与one_minus_pos是 DGP 体系中的专用原子（源码见 cvxpy/atoms/one_minus_pos.py）：diff_pos(x, y) = multiply(x, one_minus_pos(y/x))，定义域为 $x > y > 0$；one_minus_pos的定义域为 $0 < x < 1$，二者均被声明为 log-log 凹。DGP 的复合规则正是依据这些原子的 log-log 凹凸性，判定 $q$、$d$ 等复杂表达式的 log-log 凸性，从而允许它们出现在 $\leq$ 约束的左侧——这与文中第一节的数学论证完全一致。

延伸阅读与进一步探索

• 本示例出自论文Differentiating through Log-Log Convex Programs，官方文档以 derivatives 系列示例 组织，同目录下还有fundamentals.rst（可微规划基础）与structured_prediction.rst（结构化预测应用）；
• 在仓库测试中，cvxpy/tests/test_derivative.py 对backward/derivative的正确性做了大量数值校验，是理解边界行为（如不可行/无界问题的报错、非 DPP 问题的DPPError）的最佳入口；
• 关于 DGP 建模本身，doc/source/tutorial/dgp/index.rst 与 doc/source/examples/dgp/ 提供了几何规划与广义几何规划的更多案例。

实践提醒：运行本文代码前请确认已安装cvxpy及其求解器依赖（SCS 为默认支撑，requires_grad模式下 SCS 即可，无需额外安装 diffcp），并保持 numpy 可用。将 $N=2$ 扩展为更大的队列数量时，只需同步调整变量/参数的shape与对应数值即可，建模与微分逻辑完全不变。

• 科学计算