Преобразование полуопределенной программы из CVX в CVXPY со сложной целевой функцией, основанной на моделях связи

Question

Я хочу преобразовать следующий SDP (

<a href="http://oa.ee.tsinghua.edu.cn/dailinglong/publications/code/Hybrid%20precoding-based%20millimeter-wave%20massive%20MIMO-NOMA%20with%20simultaneous%20wireless%20information%20and%20power%20transfer.zip" rel="nofollow noreferrer">http://oa.ee.tsinghua.edu.cn/dailinglong/publications/code/Hybrid%20precoding-based%20millimeter-wave%20massive%20MIMO-NOMA%20with%20simultaneous%20wireless%20information%20and%20power%20transfer.zip</a>

) - который просто решает оптимизацию, подобную WMMSE - из CVX (MATLAB) в CVXPY (Python):

function [Popt, Bopt, Rsum, feasbile_point] = CVX_PA(H_eq,A,C,setf,ind,K,rf_num,sigma2,noise,P,Pconv,Pmin,Rmin)

ita = 2^Rmin-1;

cvx_begin sdp
cvx_precision best
variable powerf(K) nonnegative
variable psl(K) nonnegative
variable t1(K) nonnegative
variable t2(K) nonnegative

objective_function = myfun(powerf,t1,H_eq,A,C,setf,ind,K,rf_num,sigma2,noise);

minimize objective_function

subject to

% C2
sum(powerf) <= P;

% C3,C4
for n = 1:rf_num
    usern = setf(n,:);  % 第n个beam的用户
    usern(usern==0) = [];
    for m = 1:length(usern) 
        c3 = abs(H_eq(usern(m),usern(m)))^2*powerf(ind(n)+m)-ita*(norm(H_eq(usern(m),usern(m))))^2*sum(powerf(ind(n)+1:ind(n)+m-1));
        c4 = sigma2;
        for j = 1:rf_num               
            userj = setf(j,:);  % 第n个beam的用户
            userj(userj==0) = [];
            if j ~=  n 
                c3 = c3-ita*abs(H_eq(usern(m),userj(1)))^2*sum(powerf(ind(j)+1:ind(j+1)));
            end
            c4 = c4+abs(H_eq(usern(m),userj(1)))^2*sum(powerf(ind(j)+1:ind(j+1)));
        end
        c3 >= ita*(sigma2+noise*t1(ind(n)+m));
        c4 >= t2(ind(n)+m);
    end
end

% C5,C6,C7
for k = 1:K
    [psl(k),1;1,t1(k)] >= 0;
    [t2(k),sqrt(Pmin/Pconv);sqrt(Pmin/Pconv),1-psl(k)] >= 0;
    psl(k) <= 1;
end

cvx_end
cvx_clear

Popt = powerf;
Bopt = psl;
Rsum = -objective_function;

feasbile_point = 1;

if(strcmp(cvx_status, 'Infeasible') ||strcmp(cvx_status, 'Failed'))
    feasbile_point = 0;
end
function sume = myfun(power,t1,H_eq,A,C,setf,ind,K,rf_num,sigma2,noise)

bet = cvx(zeros(rf_num,K));
E = cvx(zeros(rf_num,K));
for n = 1:rf_num
    usern = setf(n,:);  % 第n个beam的用户
    usern(usern==0) = [];
    for m = 1:length(usern)
        bet(n,m) = (norm(H_eq(usern(m),usern(m))))^2*sum(power(ind(n)+1:ind(n)+m))+sigma2+noise*t1(ind(n)+m); % 干扰加噪声项
        for j = 1:rf_num
            if j~=n
                bet(n,m) = bet(n,m)+(norm(H_eq(usern(m),setf(j,1))))^2*sum(power(ind(j)+1:ind(j+1)));
            end
        end
        %E(n,m) = (abs(1-C(n,m)*sqrt(power(ind(n)+m))*H_eq(usern(m),usern(m))))^2+abs(C(n,m))^2*bet(n,m);
        E(n,m) = 1-2*real(C(n,m)*H_eq(usern(m),usern(m)))*sqrt(power(ind(n)+m))+abs(C(n,m))^2*bet(n,m);
    end
end

sume = 0;
for n = 1:rf_num
    usern = setf(n,:);  % 第n个beam的用户
    usern(usern==0) = [];
    for m = 1:length(usern)
        sume = sume+A(n,m)*E(n,m)/log(2)-log2(A(n,m))-1/log(2);
    end
end

Ниже мой python код

  def myfun(self,powerf,t1,H_eq, A, C, setf, ind, bet, E):
        # self, powerf, t1, H_eq, A, C, setf, ind
        self.ind = ind

        # bet = np.zeros((self.N_RF, self.K))
        # E = np.zeros((self.N_RF, self.K))
        # bet = bet1
        # E = E1
        for n in range(1, self.N_RF + 1):
            usern = setf[n - 1, :]
            usern = [item for item in usern if item != 0]
            # for m in range(1, usern.size + 1):
            # usern = np.array(usern).astype(np.int)
            usern = np.array(usern, dtype=np.int)
            for m in range(1, len(usern) + 1):
                bet[n-1, m-1] = (np.linalg.norm(H_eq[usern[m - 1]-1, usern[m - 1]-1])) ** 2 * np.sum(powerf[self.ind[n-1]:self.ind[n-1] + m - 1]) + self.sigma2 + self.noise * t1[self.ind[n-1]+m-1]
                for j in range(1, self.N_RF + 1):
                    if j != n:
                        beta[n-1, m-1] += (np.linalg.norm(H_eq[usern[m - 1], self.setf[j-1, 0]])) ** 2 * (np.sum(powerf[self.ind[j - 1]:self.ind[j]]))
                # temp1 = (np.linalg.norm(H_eq[usern[m - 1]-1, usern[m - 1]-1])) ** 2 * np.sum(powerf[self.ind[n-1]:self.ind[n-1] + m - 1]) + self.sigma2 + self.noise * t1[self.ind[n-1]+m-1]
                # for j in range(1, self.N_RF + 1):
                #     if j != n:
                #         temp1 = temp1 + (np.linalg.norm(H_eq[usern[m - 1], setf[j-1, 0]])) ** 2 * (np.sum(powerf[self.ind[j - 1]:self.ind[j]]))
                E[n-1, m-1] = 1 - 2 * (C[n-1, m-1]*H_eq[usern[m-1]-1, usern[m-1]-1]).real * np.sqrt(powerf[ind[n-1]+m-1]) + abs(C[n-1, m-1])**2 * bet[n-1, m-1]
        sume = 0
        for n in range(1, self.N_RF+1):
            usern = setf[n-1, :]
            usern = [item for item in usern if item != 0]
            for m in range(1, usern.size + 1):
                sume = sume + A[n-1, m-1] * E[n-1, m-1]/np.log(2) - np.log2(A[n-1, m-1])-1/np.log(2)
        return sume

    def CVX_PA(self, H_eq, A, C, setf, ind):
        # self, H_eq, A, C, setf, ind
        #define and solve the cvxpy problem.
        # powerf = cvx.Variable(self.K, noneg = True)
        # psl = cvx.Variable(self.K,noneg = True)
        # t1 = cvx.Variable(self.K,noneg = True)
        # t2 = cvx.Variable(self.K,noneg = True)
        powerf = cvx.Variable(self.K)
        psl = cvx.Variable(self.K)
        t1 = cvx.Variable(self.K)
        t2 = cvx.Variable(self.K)
        bet = cvx.Variable((self.N_RF, self.K))
        E = cvx.Variable((self.N_RF, self.K))
        # bet = cvx.Variable()
        # E = cvx.Variable()
        # constraints = [powerf >> 0] # only square matrix can use 只有方阵才能用>>作为正定判定
        # constraints += [psl >> 0]
        # constraints += [t1 >> 0]
        # constraints += [t2 >> 0]
        constraints = [powerf >= 0]
        constraints += [psl >= 0]
        constraints += [t1 >= 0]
        constraints += [t2 >= 0]
        # objective_func= myfun(self,powerf,t1,H_eq, A, C, setf, ind)
        objective_func = self.myfun(powerf, t1, H_eq, A, C, setf, ind,bet,E)
        for n in range(1,self.N_RF+1):
            usern = setf[n-1, :]
            usern = [item for item in userg if item != 0]
            for m in range(1,usern.size+1):
                c3 = abs(H_eq(usern[m-1],usern[m-1]))**2 * powerf[self.ind[n-1]+m-1] -self.ita * (np.linalg.norm(H_eq[usern[m-1], usern[m-1]])**2*sum(powerf[self.ind[n-1]:self.ind[n-1]+m-1]))
                c4 = sigma2
                for j in range(1,self.N_RF+1):
                    userj = setf[j-1,:]
                    userj = [item for item in userj if item !=0]
                    if j != n:
                        c3 = c3 - self.ita * abs(H_eq[usern[m-1],userj[0]])**2 * sum(powerf[self.ind[j-1]:self.ind[j]+1])
                    c4 = c4 + abs(H_eq[usern[m-1],userj[0]])**2 * sum(powerf[self.ind[j-1]:self.ind[j]+1])
                constraints += [c3 >= self.ita * (sigma2 + noise * t1[self.ind[n-1]+m-1])]
                constraints += [c4 >= t2[self.ind[n-1]+m-1]]

        constraints += [[[psl[k],1],[1,t1[k]]] >> 0 for k in range(self.K) ]
        constraints += [[[t2[k], np.sqrt(self.Pmin/self.Pconv)], [np.sqrt(self.Pmin/self.Pconv), 1-psl[k]]] >> 0 for k in range(self.K)]
        constraints += [psl[k] <= 1 for k in range[self.K]]
        constraints = [sum[popwerf] << 0,
                       ]
        prob = cvx.Problem(cvx.Minimize(objective_func), constraints)
        prob.solve()
        feas = prob.status is cvx.OPTIMAL_INACCURATE
        Popt = powerf
        Bopt = psl
        Rsum = -objective_func
        return Popt, Bopt, Rsum, feas

Однако возникают следующие ошибки:

File "*.py", line 345, in myfun
    bet[n-1, m-1] = (np.linalg.norm(H_eq[usern[m - 1]-1, usern[m - 1]-1])) ** 2 * np.sum(powerf[self.ind[n-1]:self.ind[n-1] + m - 1]) + self.sigma2 + self.noise * t1[self.ind[n-1]+m-1]
TypeError: 'Variable' object does not support item assignment