队列设计
在这个示例中，我们考虑一个 \(M/M/N\) 队列系统的设计，并研究设计对各种参数的敏感性，例如最大总延迟和总服务速率。
这个示例在论文 Differentiating through Log-Log Convex Programs 中有描述。
我们考虑 (Markovian) 队列系统的优化，其中有 \(N\) 个队列。队列系统是一组队列，其中排队的项目等待服务；排队的项目可能是操作系统中的线程，或者是网络系统的输入或输出缓冲区中的数据包。一个自然的目标是在对队列系统的各种属性施加限制的情况下，尽量减少系统的服务负载，例如最大延迟或延迟。在这个示例中，我们将这个设计问题建模为一个对数对数凸规划问题 (LLCP)，并计算设计变量对参数的敏感性。在这里考虑的队列系统被称为 \(M/M/N\) 队列。
我们假设到达第 \(i\) 个队列的项目是以速率 \(\lambda_i\) 产生的泊松过程，并且第 \(i\) 个队列的服务时间按参数 \(\mu_i\) 的指数分布进行（其中 \(i=1, \ldots, N\)）。队列系统的负载是一个函数 \(\ell: \mathbf{R}^{N}_{++} \times \mathbf{R}^{N}_{++} \to \mathbf{R}^{N}_{++}\)，其中包括到达率向量 \(\lambda\) 和服务率向量 \(\mu\)，具有以下分量:
（这是流量负载的倒数，通常用 \(\rho\) 表示。）同样，队列占用、平均延迟和系统总延迟分别是函数 \(q\)、\(w\) 和 \(d\) 关于 \(\lambda\) 和 \(\mu\) 的函数，其分量为：
其中不等式按元素逐个比较。排队系统对队列占用、平均排队延迟和总延迟设置了限制，必须满足以下条件：
其中 \(q_{\max}\), \(w_{\max}\) 和 \(d_{\max} \in \mathbf{R}^{N}_{++}\) 是参数， 不等式按元素逐个比较。此外，到达率向量 \(\lambda\) 必须至少为 \(\lambda_{\mathrm{min}} \in \mathbf{R}^{N}_{++}\)， 服务率之和不能超过 \(\mu_{\max} \in \mathbf{R}_{++}\)。
我们的设计问题是选择到达率和服务时间，以最小化服务负载的加权和 \(\gamma^T \ell(\lambda, \mu)\)， 其中 \(\gamma \in \mathbf{R}^{N}_{++}\) 是权重向量，同时满足约束条件。该问题的形式如下：
这里，\(\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 问题。 目标函数是一个多项式，约束函数 \(w\) 也是一个多项式。函数 \(d\) 和 \(q\) 都不是多项式， 但它们是对数对数凸函数；\(d\) 的对数对数凸性由组合规则得出，因为函数 \((x, y) \mapsto y - x\) 在 \(0 < x < y\) 时是对数对数凹函数，而比率函数 \((x, y) \mapsto x/y\) 是对数对数仿射函数， 在 \(y\) 上递减。通过类似的论证，可以证明 \(q\) 也是对数对数凸函数。
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')
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)
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 ]
敏感度
我们对每个变量关于参数求导。一个要点是解决方案对于 \(d_{\max}\) 和 \(\mu_{\max}\) 的值非常敏感，增加这些参数会减少服务负载，特别是对于第一个队列。
problem.solve(requires_grad=True, gp=True, eps=1e-14, max_iters=10000, mode='dense')
for var in [lam, mu, ell]:
print('变量 ', var.name())
print('对第一个分量的梯度')
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('对第二个分量的梯度')
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('')
变量 lambda
对第一个分量的梯度
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
对第二个分量的梯度
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
变量 mu
对第一个分量的梯度
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
对第二个分量的梯度
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
变量 ell
对第一个分量的梯度
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
对第二个分量的梯度
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
扰动分析
接下来，我们以一个很小的量对每个参数进行扰动，并将一阶近似预测的解与扰动问题的真实解进行比较。
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（百分比变化）：', lam_pred)
print('预测的 mu（百分比变化）：', 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（百分比变化）：', lam_actual)
print('实际的 mu（百分比变化）：', mu_actual)
预测的 lambda（百分比变化）：[2.32203282 1.77228841]
预测的 mu（百分比变化）：[1.07166961 0.94304296]
实际的 lambda（百分比变化）：[1.99504983 1.99504965]
实际的 mu（百分比变化）：[0.87148458 1.10213353]