李雅普诺夫优化方法

李雅普诺夫优化方法（Lyapunov Optimization）是一种结合稳定性分析和随机优化的框架，广泛应用于网络资源分配、动态调度和队列系统控制等领域。下面从数学角度详细介绍其核心思想和推导过程。

---

一、基本模型与李雅普诺夫函数

1. 系统动态模型

考虑一个离散时间的队列系统，系统在时隙 $t$ 的状态由队列长度向量 $\mathbf{Q}(t)=(Q_1(t),Q_2(t),\dots ,Q_N(t))$ 描述。队列动态满足：
Q i ( t + 1 ) = max ⁡ { Q i ( t ) − b i ( t ) + a i ( t ) , 0 } , ∀ i Q_i(t+1) = \max\left\{ Q_i(t) - b_i(t) + a_i(t), 0 \right\}, \quad \forall i Qi​(t+1)=max{Qi​(t)−bi​(t)+ai​(t),0},∀i
其中：

• $a_i(t)$: 时隙 $t$ 内到达队列 $i$ 的数据量，
• $b_i(t)$: 时隙 $t$ 内从队列 $i$ 服务完成的数据量。

2. 李雅普诺夫函数（Lyapunov Function）

定义二次型李雅普诺夫函数：
$$L(\mathbf{Q}(t))=\frac{1}{2}\sum _{i=1}^N{Q}_i^2(t)$$
其物理意义是衡量队列长度的“总能量”，值越大表示系统越不稳定。

---

二、李雅普诺夫漂移（Drift）

1. 单时隙漂移

定义单时隙李雅普诺夫漂移为：
$$\Delta (\mathbf{Q}(t))=\mathbb{E}\left[L(\mathbf{Q}(t+1))-L(\mathbf{Q}(t))\mid \mathbf{Q}(t)\right]$$
漂移表示队列长度变化的期望，反映系统稳定性趋势。

2. 漂移上界推导

利用队列更新方程和不等式 ( max ⁡ { x , 0 } ) 2 ≤ x 2 (\max\{x, 0\})^2 \leq x^2 (max{x,0})2≤x2，可展开：
$$\begin{array}{rl}L(\mathbf{Q}(t+1))-L(\mathbf{Q}(t))& =\frac{1}{2}\sum _{i=1}^N\left[Q_i^2(t+1)-Q_i^2(t)\right]\\ & \le \frac{1}{2}\sum _{i=1}^N\left[(Q_i(t)-b_i(t)+a_i(t){)}^2-Q_i^2(t)\right]\\ & =\frac{1}{2}\sum _{i=1}^N\left[(a_i(t)-b_i(t){)}^2+2Q_i(t)(a_i(t)-b_i(t))\right].\end{array}$$
假设到达和服务过程的二阶矩有界，即 $\mathbb{E}[(a_i(t)-b_i(t){)}^2]\le C$，则漂移满足：
$$\Delta (\mathbf{Q}(t))\le B+\sum _{i=1}^N{Q}_i(t)\mathbb{E}[a_i(t)-b_i(t)\mid \mathbf{Q}(t)],$$
其中 $B=\frac{NC}{2}$。

---

三、李雅普诺夫漂移加惩罚（Drift-Plus-Penalty）

1. 优化目标

假设系统需要最小化时间平均的代价函数 ${\lim }_{T\to \infty }\frac{1}{T}{\sum }_{t=0}^{T-1}\mathbb{E}[f(t)]$，同时保持队列稳定。构造漂移加惩罚项：
$${\Delta }_V(\mathbf{Q}(t))=\Delta (\mathbf{Q}(t))+V\cdot \mathbb{E}[f(t)\mid \mathbf{Q}(t)],$$
其中 $V>0$ 是权衡参数，调节代价与队列长度的平衡。

2. 最小化漂移加惩罚

在每个时隙 $t$，通过选择控制变量 $\mathbf{x}(t)$（如资源分配策略），求解：
$$\underset{\mathbf{x}(t)\in {\mathcal{X}}_t}{\min }\left\{\sum _{i=1}^N{Q}_i(t)\mathbb{E}[a_i(t)-b_i(t)\mid \mathbf{Q}(t)]+V\cdot \mathbb{E}[f(t)\mid \mathbf{Q}(t)]\right\},$$
其中 ${\mathcal{X}}_t$ 是可行控制集合。

---

四、性能分析

1. 稳定性保证

如果漂移加惩罚满足：
$${\Delta }_V(\mathbf{Q}(t))\le B-ϵ\sum _{i=1}^N{Q}_i(t)+V{f}^{\ast },$$
其中 $ϵ>0$，$f^{\ast }$ 为理论最优代价，则可证明队列均值为 $O(V)$，即：
$$\underset{T\to \infty }{\mathrm{limsup}}\frac{1}{T}\sum _{t=0}^{T-1}\sum _{i=1}^N\mathbb{E}[Q_i(t)]\le \frac{B+V(f^{\ast }-\bar{f})}{ϵ},$$
其中 $\bar{f}$ 为实际代价。

2. 最优性间隙

时间平均代价满足：
$$\underset{T\to \infty }{\mathrm{limsup}}\frac{1}{T}\sum _{t=0}^{T-1}\mathbb{E}[f(t)]\le f^{\ast }+\frac{B}{V}.$$
增大 $V$ 可逼近最优解，但队列长度随 $V$ 线性增长。

---

五、实例：动态资源分配

问题：最小化平均功耗，同时稳定处理队列。

1. 模型设定：

   • 控制变量：功率分配 $p(t)\in [0,P_{\max }]$,
   • 服务率：$b(t)=\log (1+p(t)h(t))$,
   • 到达率：$a(t)\sim \text{Poisson}(\lambda )$,
   • 代价函数：$f(t)=p(t)$.

2. 策略设计：
   每个时隙求解：
   min ⁡ p ( t ) { Q ( t ) E [ a ( t ) − log ⁡ ( 1 + p ( t ) h ( t ) ) ] + V p ( t ) } . \min_{p(t)} \left\{ Q(t) \mathbb{E}[a(t) - \log(1 + p(t) h(t))] + V p(t) \right\}. p(t)min​{Q(t)E[a(t)−log(1+p(t)h(t))]+Vp(t)}.
   前面一项是为了稳定队列，后面一项是最小化代价。

---

六、总结

李雅普诺夫优化的核心步骤：

1. 定义李雅普诺夫函数衡量系统稳定性。
2. 推导漂移上界，结合惩罚项构造优化目标。
3. 在每个时隙最小化漂移加惩罚项，得到在线策略。
4. 通过参数 $V$ 调节稳定性和最优性的权衡。

优势：

• 将复杂随机优化问题分解为单时隙优化。
• 不依赖先验统计信息，适用于非稳态系统。
• 理论保证稳定性和最优性间隙。

公式总结：

• 队列更新： Q i ( t + 1 ) = max ⁡ { Q i ( t ) − b i ( t ) + a i ( t ) , 0 } Q_i(t+1) = \max\{Q_i(t) - b_i(t) + a_i(t), 0\} Qi​(t+1)=max{Qi​(t)−bi​(t)+ai​(t),0}
• 李雅普诺夫函数：$L(\mathbf{Q}(t))=\frac{1}{2}\sum Q_i^2(t)$
• 漂移加惩罚：${\Delta }_V(\mathbf{Q}(t))=\Delta (\mathbf{Q}(t))+V\cdot \mathbb{E}[f(t)]$

通过上述框架，李雅普诺夫优化为动态随机系统提供了一种优雅的联合稳定性和性能优化的方法。

---

多用户动态资源分配的数学建模

---

1. 系统模型

考虑一个蜂窝基站为 $N$ 个用户提供服务的离散时隙系统，每个时隙 t ∈ { 0 , 1 , 2 , . . . , T } t \in \{0,1,2,...,T\} t∈{0,1,2,...,T} 包含以下动态过程：

1. 信道状态：用户 $i$ 的时变信道增益为 $h_i(t)\sim \text{Rayleigh}(h_{\text{mean},i})$
2. 数据到达：用户 $i$ 的数据包到达过程服从泊松分布 $a_i(t)\sim \text{Poisson}({\lambda }_i)$
3. 队列动态：用户 $i$ 的队列长度 $Q_i(t)$ 更新方程为：
   Q i ( t + 1 ) = max ⁡ { Q i ( t ) − b i ( t ) + a i ( t ) ,   0 } Q_i(t+1) = \max\left\{ Q_i(t) - b_i(t) + a_i(t),\ 0 \right\} Qi​(t+1)=max{Qi​(t)−bi​(t)+ai​(t), 0}
   其中 $b_i(t)$ 为服务速率，由功率分配决定。
4. 功率约束：基站总发射功率满足 ${\sum }_{i=1}^N{p}_i(t)\le P_{\max }$，且 $p_i(t)\ge 0$。

---

2. 控制目标

设计动态功率分配策略 { p i ( t ) } \{p_i(t)\} {pi​(t)} 以：

1. 维持队列稳定：${\lim }_{T\to \infty }\frac{1}{T}{\sum }_{t=0}^{T-1}\mathbb{E}[Q_i(t)]<\infty ,\forall i$
2. 最小化长期平均功率：$\min {\lim }_{T\to \infty }\frac{1}{T}{\sum }_{t=0}^{T-1}\mathbb{E}\left[{\sum }_{i=1}^N{p}_i(t)\right]$

---

3. 李雅普诺夫优化框架

3.1 李雅普诺夫函数

定义二次李雅普诺夫函数表征队列积压：
$$L(\mathbf{Q}(t))=\frac{1}{2}\sum _{i=1}^N{Q}_i(t{)}^2$$

3.2 条件李雅普诺夫漂移

定义单时隙漂移量：
$$\Delta (\mathbf{Q}(t))=\mathbb{E}\left[L(\mathbf{Q}(t+1))-L(\mathbf{Q}(t))\mid \mathbf{Q}(t)\right]$$

3.3 漂移加惩罚项

引入权衡参数 $V>0$，将功率优化融入目标：
$$\Delta (\mathbf{Q}(t))+V\cdot \mathbb{E}\left[\sum _{i=1}^N{p}_i(t)\mid \mathbf{Q}(t)\right]$$

---

4. 问题转化

4.1 漂移上界推导

利用队列更新方程展开 $L(\mathbf{Q}(t+1))$，可得不等式：
$$\Delta (\mathbf{Q}(t))\le B+\sum _{i=1}^N{Q}_i(t)\mathbb{E}[a_i(t)-b_i(t)\mid \mathbf{Q}(t)]$$
其中 $B=\frac{1}{2}{\sum }_{i=1}^N\left(a_i^2+b_i^2\right)$ 为有界常量项。

4.2 优化问题设计

每个时隙 $t$ 求解以下问题以最小化漂移加惩罚项：
$$\underset{p_i(\cdot )}{\min }\sum _{i=1}^N\left(-Q_i(t)b_i(t)+V{p}_i(t)\right)$$
约束条件：
$$\{\begin{array}{l}{\sum }_{i=1}^N{p}_i(t)\le P_{\max }\\ p_i(t)\ge 0,\forall i\\ b_i(t)={\log }_2\left(1+\frac{p_i(t)h_i(t)}{N_0}\right)\text{(香农容量)}\end{array}$$

---

5. 凸优化形式

将服务速率 $b_i(t)$ 表达式代入目标函数，问题等价于：
min ⁡ { p i } ∑ i = 1 N ( − Q i ( t ) log ⁡ 2 ( 1 + p i h i ( t ) N 0 ) + V p i ) \min_{\{p_i\}} \sum_{i=1}^N \left( -Q_i(t) \log_2\left(1 + \frac{p_i h_i(t)}{N_0}\right) + V p_i \right) {pi​}min​i=1∑N​(−Qi​(t)log2​(1+N0​pi​hi​(t)​)+Vpi​)
Subject to:
$$\{\begin{array}{l}{\sum }_{i=1}^N{p}_i\le P_{\max }\\ p_i\ge 0,\forall i\end{array}$$

---

6. 模型特性

1. 凸性保证：目标函数为凸函数（$\log (1+x)$ 的凹性保留在求和后的线性组合中），约束为线性，属凸优化问题
2. 实时性：仅依赖当前时隙的队列状态 $Q_i(t)$ 与信道状态 $h_i(t)$
3. 稳定性证明：通过李雅普诺夫方法可证，当到达率在系统容量域内时，队列稳定且平均功率收敛到最优值的邻域内

---

7. 关键公式总结

组件	数学表达式
队列更新	$ Q_i(t+1) = \max\left{ Q_i(t) - b_i(t) + a_i(t), 0 \right} $
服务速率	$ b_i(t) = \log_2\left(1 + \frac{p_i(t) h_i(t)}{N_0}\right) $
优化目标	$ \min \sum_{i=1}^N \left( -Q_i(t) b_i(t) + V p_i(t) \right) $
功率约束	$ \sum_{i=1}^N p_i(t) \leq P_{\max} $

---

此模型通过动态权衡队列积压与功率消耗，实现了在未知随机信道和随机到达下的最优资源分配。代码实现严格遵循上述数学框架，通过凸优化求解器在每个时隙实时计算最优功率分配。

import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp

# ================= 系统参数 ======================
N_users = 3  # 用户数量
T = 500  # 时隙总数
V = 15  # 李雅普诺夫权衡参数
P_max_total = 5  # 总功率约束 (W)
noise_power = 0.1  # 噪声功率 (W)
arrival_rates = [2, 3, 1.5]  # 各用户到达率
h_means = [0.8, 1.2, 1.0]  # 信道均值

# ================= 状态初始化 ====================
Qs = np.zeros(N_users)
power_history = np.zeros((T, N_users))
queue_history = np.zeros((T, N_users))
h_history = np.zeros((T, N_users))
b_history = np.zeros((T, N_users))  # 新增传输数据量记录


# ================= 辅助函数 ======================
def multi_channel():
    return [np.abs(np.random.normal(loc=mean, scale=0.3 * mean)) for mean in h_means]


def multi_arrival():
    return [np.random.poisson(lam) for lam in arrival_rates]


def service_rates(powers, h):
    return [np.log2(1 + p * h[i] / noise_power) for i, p in enumerate(powers)]


# ================= 核心优化器 ====================
def lyapunov_optimization_multi(Qs, h, V):
    p = cp.Variable(N_users)
    obj_terms = []
    for i in range(N_users):
        snr = p[i] * h[i] / noise_power
        rate = cp.log(1 + snr) / np.log(2)
        obj_terms.append(Qs[i] * (-rate) + V * p[i])

    constraints = [p >= 0, cp.sum(p) <= P_max_total]

    prob = cp.Problem(cp.Minimize(sum(obj_terms)), constraints)
    prob.solve(solver=cp.ECOS)

    return np.clip(p.value, 0, P_max_total) if prob.status == 'optimal' else np.zeros(N_users)


# ================= 主循环 ========================
for t in range(T):
    h = multi_channel()
    a_t = multi_arrival()
    h_history[t] = h

    p_alloc = lyapunov_optimization_multi(Qs, h, V)
    b_t = service_rates(p_alloc, h)

    Qs = np.maximum(Qs - np.array(b_t) + np.array(a_t), 0)

    # 记录所有状态
    power_history[t] = p_alloc
    queue_history[t] = Qs
    b_history[t] = b_t  # 新增记录

# ================= 可视化增强 ======================
plt.figure(figsize=(16, 12))

# 功率分配三维图
ax1 = plt.subplot(421, projection='3d')
for i in range(N_users):
    ax1.plot(range(T), [i] * T, power_history[:, i], label=f'User {i + 1}')
ax1.set_title('3D Power Allocation')
ax1.set_zlabel('Power (W)')

# 队列堆积图
plt.subplot(422)
plt.stackplot(range(T), queue_history.T, labels=[f'User {i + 1}' for i in range(N_users)])
plt.title('Cumulative Queue Backlog')
plt.ylabel('Packets')

# 传输数据量可视化（新增）
plt.subplot(423)
for i in range(N_users):
    plt.plot(b_history[:, i], alpha=0.8, linestyle='--', linewidth=1, label=f'User {i + 1}')
plt.title('Service Rate (Transmitted Data)')
plt.ylabel('Bits/slot')
plt.legend()
plt.grid(True)

# 总功率监控
plt.subplot(424)
plt.plot(power_history.sum(axis=1), color='darkred')
plt.axhline(P_max_total, ls='--', color='grey', label='Power Limit')
plt.title('Total Power Consumption')
plt.xlabel('Time Slot')
plt.legend()

# 信道质量热力图（新增辅助分析）
plt.subplot(425)
plt.imshow(h_history.T, aspect='auto', cmap='viridis')
plt.colorbar(label='Channel Gain')
plt.title('Channel Quality Heatmap')
plt.ylabel('User Index')
plt.xlabel('Time Slot')

# 到达与服务量对比（新增）
plt.subplot(426)
for i in range(N_users):
    plt.plot(np.cumsum([a[i] for a in np.array([multi_arrival() for _ in range(T)])]),
             label=f'Arrival {i + 1}', linestyle=':')
    plt.plot(np.cumsum(b_history[:, i]),
             label=f'Service {i + 1}', linewidth=1.5)
plt.title('Cumulative Arrival vs Service')
plt.legend()

plt.tight_layout()
plt.show()