How to implement the Linear Relaxation (LR) algorithm?

[zlkmn]

How to implement the Linear Relaxation (LR) algorithm? After the decision variable x is relaxed to [0,1], problem P1 does not conform to the DCP rule in CVXPY (Disciplined Convex Programming), and cannot be solved by CVXPY package.

import numpy as np
import cvxpy as cp
import scipy.io as sio  # import scipy.io for .mat file I/
import time

def lr_algorithm(h, weights=[]):
    """
    线性松弛算法实现（基于CVXPY）
    :param h: 信道增益数组
    :param weights: 设备权重列表
    :return: (最大计算速率, 最优卸载决策, 能量传输时间, 传输时间分配)
    """
    N = len(h)  # 设备数量

    # ===== 系统参数 =====
    mu = 0.7  # 能量收集效率
    P = 3  # AP发射功率 (W)
    B = 2e6  # 带宽 (Hz)
    N0 = 1e-10  # 噪声功率 (W)
    k = 1e-26  # 本地计算能耗系数
    phi = 100  # 计算周期数/bit
    vu = 1.1  # 数据上传开销系数

    if len(weights) == 0:  # 默认权重
        weights = [1.5 if i % 2 else 1 for i in range(N)]
    weights = np.array(weights)

    # ===== 定义优化变量 =====
    x = cp.Variable(N, nonneg=True)  # 松弛后的卸载决策 [0,1]
    tau = cp.Variable(N, nonneg=True)  # 传输时间分配
    a = cp.Variable(nonneg=True)  # 能量传输时间

    # ===== 目标函数：加权计算速率和 =====
    eta1 = (mu * P) ** (1 / 3) / phi
    objective = 0

    # 显式遍历每个设备，避免向量化潜在问题
    for i in range(N):
        # 本地计算速率
        r_L = eta1 * (h[i] / k) ** (1 / 3) * a ** (1 / 3)
        # 云端计算速率
        SNR = (mu * P * h[i] ** 2 * a) / N0
        r_O = (B * tau[i] / vu) * cp.log(1 + SNR / (tau[i] + 1e-10)) / np.log(2)
        # 加权累加
        objective += weights[i] * ((1 - x[i]) * r_L + x[i] * r_O)

    # ===== 约束条件 =====
    constraints = [
        a + cp.sum(tau) <= 1,  # 总时间约束
        a >= 0,
        tau >= 0,
        x <= 1,  # 松弛变量上界
        tau <= x  # 传输时间约束（x_i=0时tau_i=0）
    ]

    # ===== 求解凸问题 =====
    problem = cp.Problem(cp.Maximize(objective), constraints)
    print("r_L DCP:", r_L.is_dcp())
    print("SNR:", SNR.is_dcp())
    print("r_O DCP:", r_O.is_dcp())
    print("Objective DCP:", objective.is_dcp())
    problem.solve(solver=cp.ECOS, verbose=False)

    # ===== 二值化处理 =====
    x_binary = np.zeros(N)
    if a.value is None:  # 处理求解失败的情况
        return 0, x_binary, 0, np.zeros(N)

    for i in range(N):
        # 计算本地和云端速率
        r_L_val = eta1 * (h[i] / k) ** (1 / 3) * a.value ** (1 / 3)
        if tau[i].value is None or tau[i].value <= 1e-10:
            r_O_val = 0
        else:
            SNR_val = (mu * P * h[i] ** 2 * a.value) / (tau[i].value + 1e-10)
            r_O_val = (B * tau[i].value / vu) * np.log2(1 + SNR_val / N0)
        x_binary[i] = 1 if r_O_val >= r_L_val else 0

    # ===== 重新求解资源分配 =====
    active_devices = np.where(x_binary == 1)[0]
    if len(active_devices) == 0:
        return 0, x_binary, a.value, np.zeros(N)

    a_p = cp.Variable(nonneg=True)
    tau_p = cp.Variable(len(active_devices), nonneg=True)

    obj_act = 0
    # 本地设备贡献
    for i in range(N):
        if x_binary[i] == 0:
            r_L = eta1 * (h[i] / k) ** (1 / 3) * a_p ** (1 / 3)
            obj_act += weights[i] * r_L
    # 卸载设备贡献
    for idx, j in enumerate(active_devices):
        SNR = (mu * P * h[j] ** 2 * a_p) / (tau_p[idx] + 1e-10)
        r_O = (B * tau_p[idx] / vu) * cp.log(1 + SNR / N0) / np.log(2)
        obj_act += weights[j] * r_O

    constraints_act = [
        a_p + cp.sum(tau_p) <= 1,
        a_p >= 0,
        tau_p >= 0
    ]

    problem_act = cp.Problem(cp.Maximize(obj_act), constraints_act)
    problem_act.solve(solver=cp.ECOS, verbose=False)

    # ===== 组装最终结果 =====
    final_tau = np.zeros(N)
    if tau_p.value is not None:
        final_tau[active_devices] = tau_p.value
    final_rate = problem_act.value if problem_act.status == 'optimal' else 0

    return final_rate, x_binary, a_p.value if a_p.value is not None else 0, final_tau

# 测试代码
if __name__ == "__main__":
    # test all data
    # ===== 系统参数设置 =====

    K = [10]  # number of users
    N = 6000  # number of channel # 每个用户数量下的测试样本数（信道实例数）


    for k in K:
        # Load data
        channel = sio.loadmat('./data/data_%d' % int(k))['input_h']
        gain = sio.loadmat('./data/data_%d' % int(k))['output_obj']

        start_time = time.time()
        gain_his = []  # 存储每个信道实例的计算速率
        gain_his_ratio = []  # 存储归一化后的速率（当前速率 / 最优速率）
        mode_his = []  # 存储最优卸载决策
        for i in range(N):
            if i % (N // 10) == 0:
                print("%0.1f" % (i / N))  # 打印进度（每处理10%输出一次）

            i_idx = i+24000

            h = channel[i_idx, :]  # 获取第i个信道实例的增益

            # the LR method
            rate, mode, a, tau = lr_algorithm(h)

            # memorize the largest reward
            gain_his.append(rate)
            gain_his_ratio.append(gain_his[-1] / gain[i_idx][0])

            mode_his.append(mode)

        total_time = time.time() - start_time
        print('time_cost:%s' % total_time)
        print('average time per channel:%s' % (total_time / N))
        print("Max computation rate:", max(gain_his))
        print("Average computation rate:", sum(gain_his)/N)
        print("gain/max ratio: ", sum(gain_his_ratio) / N)

Error：

r_L` DCP: True
SNR: True
r_O DCP: False
Objective DCP: False
    raise DCPError(
cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP.