CVXPY DCP规则不满足报错求助：无人机通信优化算法复现遇阻

问题背景 #

复现论文《Energy Minimization for Wireless Communication With Rotary-Wing UAV》的无人机通信优化算法，采用CVXPY库结合SCA算法，目标为最小化无人机能耗的同时最大化通信速率，输入数据来自coordinates.csv文件。尝试用泰勒近似处理非凸函数后，仍出现"DCP Problem does not follow DCP rules"报错，相关代码如下：

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

# 读取坐标文件
df = pd.read_csv('coordinates.csv')
qk = df.to_numpy()
qk = qk[1:]
wk = qk.copy()

# 参数设置
Ph = 168.5
Pc = 5
epsilon = 1e-6  # 收敛阈值
max_dist = 105.57
referenceSnr = 177827.94
alpha = 1
datasize = 10
E0 = 36.76
max_iter = 10

# 初始化变量
qk_tilt = cp.Variable((9,2), nonneg=True)
# 错误：将迭代参考点定义为CVXPY变量，而非固定常数
qk_tilt_l = cp.Variable((9, 2), value=(wk))
nk = cp.Variable(9, nonneg=True)
x = cp.Variable(8)
index = 0

# SCA迭代
while index <1:
    constraints = [cp.sum(x) <= max_dist**2]
    constraints += [cp.sum_squares(qk_tilt[i+1]-qk_tilt[i]) <= x[i] for i in range(8)]
    
    # 计算当前参考点的速率值
    Rk = [cp.log(1+((referenceSnr)/(100**2+cp.sum_squares(qk_tilt_l[i]-wk[i]))))/cp.log(2) for i in range(9)]
    # 泰勒近似系数计算
    bk = [((cp.log(cp.exp(1))/cp.log(2)) * referenceSnr * alpha) / (100 ** 2 + cp.sum_squares(qk_tilt_l[i] - wk[i]) * ((100**2 + cp.sum_squares(qk_tilt_l[i] - wk[i])) + referenceSnr)) for i in range(9)]
    term1 = Rk.copy()
    term2 = [bk[i] * (cp.sum_squares(qk_tilt[i]-wk[i])-cp.sum_squares(qk_tilt_l[i]-wk[i])) for i in range(9)]
    
    constraints += [nk[i] <=(term1[i] - term2[i]) for i in range(9)]
    constraints += [nk[i] >= 0 for i in range(9)]
    
    # 冗余处理：nk已为非负变量，无需cp.pos
    nk_arr = cp.vstack([cp.pos(nk[i]) for i in range(9)])
    obj = cp.Minimize(E0*cp.sum(x)+cp.sum((Ph+Pc)*datasize*nk_arr))
    
    prob = cp.Problem(obj, constraints)
    try:
        prob.solve(verbose=True)
        if cp.sum_squares(qk_tilt_l - qk_tilt) < epsilon**2:
            break
    except Exception as e :
        print("ERROR",e) 
    
    qk_tilt_l.value = qk_tilt.value
    l = l + 1 
    index+=1

# 绘图
plt.figure()
plt.scatter(wk[:, 0], wk[:, 1], label='wk')
plt.scatter(qk_tilt_l.value[:, 0], qk_tilt_l.value[:, 1], label='qk_tilt_l')
plt.legend()
plt.show()

报错原因分析 #

1. 迭代参考点变量类型错误：qk_tilt_l是SCA迭代中每一轮的固定参考点，应定义为numpy常数数组，但代码中将其设为CVXPY优化变量，导致约束中出现变量间的非线性组合，违反DCP规则。
2. 泰勒近似符号与凸性不匹配：原速率函数是关于无人机位置的凹函数，泰勒近似的方向错误，导致约束项不满足DCP的凸/凹性要求。
3. 冗余变量处理：nk已声明为非负变量，使用cp.pos(nk[i])属于冗余操作，且cp.vstack会导致目标函数的表达式结构不符合DCP规范。
4. 约束缺失：x表示相邻无人机位置的距离平方，未添加x[i] >= 0的非负约束，可能引发数值问题。
5. 迭代逻辑错误：l变量未初始化就进行自增操作，且循环条件index <1仅执行一次迭代，不符合SCA多轮迭代的要求。

解决方案 #

1. 修正迭代参考点类型：将qk_tilt_l改为numpy数组，存储上一轮迭代的最优位置，作为当前轮的固定常数。
2. 修正泰勒近似表达式：根据凹函数的一阶泰勒展开性质，调整约束项的符号，确保约束符合DCP规则。
3. 简化目标函数：直接对nk求和，去掉冗余的cp.pos和cp.vstack操作。
4. 补充约束：添加x[i] >= 0的非负约束。
5. 修复迭代逻辑：初始化迭代计数器，设置合理的迭代终止条件（达到最大迭代次数或满足收敛阈值）。

修正后代码 #

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

# 读取坐标文件
df = pd.read_csv('coordinates.csv')
qk = df.to_numpy()
qk = qk[1:]
wk = qk.copy()
num_uav = wk.shape[0]

# 参数设置
Ph = 168.5
Pc = 5
epsilon = 1e-6  # 收敛阈值
max_dist = 105.57
referenceSnr = 177827.94
alpha = 1
datasize = 10
E0 = 36.76
max_iter = 10

# 初始化变量
qk_tilt = cp.Variable((num_uav, 2), nonneg=True)
# 修正：用numpy数组存储迭代参考点，初始值为wk
qk_tilt_l = wk.copy().astype(np.float64)
nk = cp.Variable(num_uav, nonneg=True)
x = cp.Variable(num_uav - 1, nonneg=True)  # 添加非负约束

# SCA迭代
l = 0
converged = False
while l < max_iter and not converged:
    constraints = []
    # 总飞行距离平方约束
    constraints.append(cp.sum(x) <= max_dist ** 2)
    # 相邻UAV距离平方约束
    for i in range(num_uav - 1):
        constraints.append(cp.sum_squares(qk_tilt[i+1] - qk_tilt[i]) <= x[i])
    
    # 计算当前参考点的速率值和泰勒近似系数
    Rk = []
    grad_Rk = []
    for i in range(num_uav):
        dist_sq = np.sum((qk_tilt_l[i] - wk[i]) ** 2)
        denom = 100**2 + dist_sq
        # 当前参考点的速率值（常数）
        r_val = np.log(1 + referenceSnr / denom) / np.log(2)
        Rk.append(r_val)
        # 计算速率函数对qk_tilt[i]的梯度（常数向量）
        grad = (np.log(np.e)/np.log(2)) * referenceSnr * (-2)*(qk_tilt_l[i] - wk[i]) / (denom * (denom + referenceSnr))
        grad_Rk.append(grad)
    
    # 构建SCA约束：nk[i] <= Rk(q^t) + ∇Rk(q^t)^T(q - q^t)
    # 因为Rk是凹函数，一阶泰勒是上界，约束nk <= 上界等价于用凸近似替代原凹约束
    for i in range(num_uav):
        linear_term = grad_Rk[i].T @ (qk_tilt[i] - qk_tilt_l[i])
        constraints.append(nk[i] <= Rk[i] + linear_term)
    
    # 定义目标函数：最小化能耗
    obj = cp.Minimize(E0 * cp.sum(x) + (Ph + Pc) * datasize * cp.sum(nk))
    prob = cp.Problem(obj, constraints)
    
    try:
        prob.solve(verbose=True)
        # 检查收敛：位置变化的平方和小于阈值
        if np.sum((qk_tilt.value - qk_tilt_l) ** 2) < epsilon ** 2:
            converged = True
        # 更新参考点
        qk_tilt_l = qk_tilt.value.copy()
        l += 1
    except Exception as e:
        print("ERROR:", e)
        break

# 绘图
plt.figure()
plt.scatter(wk[:, 0], wk[:, 1], label='初始位置wk')
plt.scatter(qk_tilt_l[:, 0], qk_tilt_l[:, 1], label='优化后位置qk_tilt')
plt.xlabel('X坐标')
plt.ylabel('Y坐标')
plt.legend()
plt.grid(True)
plt.show()

内容的提问来源于stack exchange，提问作者khebila zaurt