科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载导读本文基于 CVXPY 官方示例 queuing_design.rst完整讲解如何将 M/M/N 排队系统的服务负载最小化问题建模为 log-log 凸规划LLCP并使用 CVXPY 的requires_gradTrue可微求解能力通过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$ 的指数分布$i1,\ldots,N$。系统核心度量是服务负载函数 $\ell: \mathbf{R}^{N}{} \times \mathbf{R}^{N}{} \to \mathbf{R}^{N}_{}$$$ \ell_i(\lambda, \mu) \frac{\mu_i}{\lambda_i}, \quad i1, \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_{i1}^{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.pydiff_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(posTrue, shape(2,), namemu) lam cp.Variable(posTrue, shape(2,), namelambda) ell cp.Variable(posTrue, shape(2,), nameell) w_max cp.Parameter(posTrue, shape(2,), valuenp.array([2.5, 3.0]), namew_max) d_max cp.Parameter(posTrue, shape(2,), valuenp.array([2., 2.]), named_max) q_max cp.Parameter(posTrue, shape(2,), valuenp.array([4., 5.0]), nameq_max) lam_min cp.Parameter(posTrue, shape(2,), valuenp.array([0.5, 0.8]), namelambda_min) mu_max cp.Parameter(posTrue, value3.0, namemu_max) gamma cp.Parameter(posTrue, shape(2,), valuenp.array([1.0, 2.0]), namegamma)这里的关键点所有变量均以posTrue声明强制取正对应数学定义域 $\mathbf{R}^{N}_{}$所有参数也以posTrue声明同样约束为正向量参数使用shape(2,)本示例取 $N2$ 个队列标量参数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_gradTrue, gpTrue, eps1e-14, max_iters10000, modedense)求解结果目标函数最优值为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_gradTrue, gpTrue, ...)的求解器调用细节gpTrue表示按 DGP/LLCP 模式解析问题requires_gradTrue则要求在求解时缓存可微求解器的内部数据供后续backward()/derivative()使用。从 cvxpy/problems/problem.py 的源码可以看到启用requires_grad的三条硬性约束问题必须是DPP 合规的problem.is_dpp(dpp_context)dpp_context在gpTrue时为dgp否则抛出DPPError求解器必须为SCS 或 DIFFCP否则抛出ValueError若选用DIFFCP需要安装独立的diffcpPython 包。eps1e-14与max_iters10000是透传给 SCS 的收敛容差与迭代上限modedense指定使用稠密后端。从实现上看requires_gradTrue的求解结果会被缓存到self._solver_cache[s.DIFFCP]中backward_cache其中保存了求解器导出的微分算子D/DT这正是后续两类微分 API 的数据基础。灵敏度分析用backward()计算变量对参数的梯度在 LLCP 设计问题中我们关心的核心问题是最优设计对哪些参数最敏感例如放宽最大总延迟 $d_{\max}$ 或提高服务率总预算 $\mu_{\max}$ 对服务负载分别有多大影响这可以通过backward()求解。problem.solve(requires_gradTrue, gpTrue, eps1e-14, max_iters10000, modedense) 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_gradTrue, ...)之后调用且要求问题状态为有解否则抛出SolverError。扰动分析用derivative()做一阶预测并验证灵敏度分析回答的是参数变化一点点解怎么变的线性化问题扰动分析则进一步用一阶近似预测参数扰动后的新解并与重新求解的真实解对比验证线性近似的精度。problem.solve(requires_gradTrue, gpTrue, eps1e-14, max_iters10000, modedense) 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, gpTrue, eps1e-14, max_iters10000) 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]这里的操作协议是给每个参数设置param.delta param.value * delta即相对扰动 $1%$调用problem.derivative()它会填充变量delta属性——即一阶预测的变量变化量把扰动真正施加到参数上param.value param.delta重新求解原始问题得到真实解分别计算预测与实际的百分比变化并对比。对比可见一阶预测与真实变化的量级完全一致约 $1%\sim 2%$ 量级方向相同、数值接近验证了灵敏度分析结论的可靠性当然在 $1%$ 的有限扰动下也存在可观测的非线性偏差这正是线性近似的固有误差。源码级机制problem.derivative() 是backward()的正向模式对应物——它把参数扰动 $\delta p$ 通过解映射的导数映射为变量最优值的扰动 $\delta x$。实现要点未设置delta的参数默认扰动为 0np.broadcast_to(0.0, param.shape)沿约简链正向调用param_forward把参数扰动传导到参数化问题利用缓存的D算子计算内部解扰动再经逆序的var_forward还原为外层变量delta副作用是填充每个Variable的delta属性同样要求先以requires_gradTrue求解且问题有解。这两类 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_gradTrue时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.pydiff_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 可用。将 $N2$ 扩展为更大的队列数量时只需同步调整变量/参数的shape与对应数值即可建模与微分逻辑完全不变。赞分享科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载相关推荐algo 项目实战归并排序与快速排序的 O(n log n) 实现、复杂度分析与第 K 大元素查找algo 项目实战归并排序与快速排序的 O n log n 实现、复杂度分析与第 K 大元素查找 本篇文章以仓库笔记 notes/12_sorts/readm示例工程LeetCode 0004 寻找两个正序数组的中位数从暴力合并到 O(log(min(n,m))) 二分划分含多语言实现LeetCode 0004 寻找两个正序数组的中位数从暴力合并到 O log min n,m 二分划分含多语言实现 导读 本篇围绕 LeetCode 第示例工程教程30 Seconds of Code 中的归并排序分治思想、递归实现与 O(n log n) 复杂度剖析30 Seconds of Code 中的归并排序分治思想、递归实现与 O n log n 复杂度剖析 归并排序Merge Sort是一种高效的、基于比较教程文档上一篇C-Shopping组件设计模式可复用UI组件库的构建思路下一篇如何在Unity中实现智能AIGoAP目标导向行为规划完整教程创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考