function [S_m_out, fval_2, flag_2] = S_m_calc_CVX(alpha, S_m_r, p_m, q_m_x, q_m_y)
%%参数定义
global P_charge;
global H;
global d_min;
global V_max;
global S;
global P_max;
global freq; %载波频率
global light_speed; %光速
global B;
global yeta_los;
global delta2;
global K;
global M;
global N;
global R_th;
global P_user;
%%变量初始化
A_m_k = zeros(N,M,K);
for m=1:1:M %%计算 c_m_k'(n)
for k=1:1:K
for n=1:1:N
dist2 = H*H + (q_m_x(n,m)-P_user(1,k))^2 + (q_m_y(n,m)-P_user(1,k))^2;
L_current = (4*pi*freq/light_speed)^2 * dist2 * yeta_los;
P_receive_current = p_m(n,m)/L_current;
P_receive_total = 0;
for m1=1:1:M;
dist2 = H*H + (q_m_x(n,m1)-P_user(1,k))^2 + (q_m_y(n,m1)-P_user(1,k))^2;
L_current = (4*pi*freq/light_speed)^2 * dist2 * yeta_los;
P_receive_temp = p_m(n,m1)/L_current;
P_receive_total = P_receive_total + P_receive_temp;
end
A_m_k(n,m,k) = P_receive_current/(P_receive_total - P_receive_current + B*delta2);
end
end
end
alpha_1 = reshape(alpha, N*M, K);
alpha_m_n_0 = sum(alpha_1,2);
A_m_k_1 = reshape(A_m_k, N*M, K);
S_m_r_0 = S_m_r;
for k=1:1:K-1
S_m_r_0 = cat(3,S_m_r_0,S_m_r);
end
alpha_m_n = alpha_m_n_0;
for k=1:1:K-1
alpha_m_n = cat(2,alpha_m_n,alpha_m_n_0);
end
S_m_r_1 = reshape(S_m_r_0, N*M, K);
% CC = -log(S_m_r_1)./log(2) + 1/log(2);
CC = log(S_m_r_1)./log(2);
cvx_begin
cvx_solver mosek
variable S_m_v(N,M)
expression S_m_v_cat
S_m_v_cat = S_m_v;
for k=1:1:K-1 %扩展S_m_r的维度到N*M*K,维度统一方便后续CVX计算
S_m_v_cat = cat(3,S_m_v_cat,S_m_v);
end
S_m_v_cat = reshape(S_m_v_cat, N*M, K);
maximize(sum(sum(alpha_1.*(-rel_entr(1, S_m_v_cat+A_m_k_1)./log(2) - (CC + (S_m_v_cat - S_m_r_1) ./(S_m_r_1.*log(2))) ))).*B./(M*N));
subject to
for k=1:1:K
sum(alpha_1(:,k).*(-rel_entr(1, S_m_v_cat(:,k)+A_m_k_1(:,k))./log(2) - (CC(:,k) + (S_m_v_cat(:,k) - S_m_r_1(:,k)) ./(S_m_r_1(:,k).*log(2))) )).*B./N >= R_th;
end
S_m_v >= sum(alpha,3);
S_m_v <= S;
cvx_end
S_m_out = reshape(S_m_v(:,1),N,M);
fval_2 = cvx_optval;
flag_2 = cvx_status;
The output shows:
CVX Warning:
Models involving "rel_entr" or other functions in the log, exp, and entropy
family are solved using an experimental successive approximation method.
This method is slower and less reliable than the method CVX employs for
other models. Please see the section of the user's guide entitled
The successive approximation method
for more details about the approach, and for instructions on how to
suppress this warning message in the future.
=====================================
Using Pade approximation for exponential
cone with parameters m=3, k=3
=====================================
=====================================
Using Pade approximation for exponential
cone with parameters m=3, k=3
=====================================
=====================================
Using Pade approximation for exponential
cone with parameters m=3, k=3
=====================================
Calling Mosek 8.0.0.60: 105 variables, 22 equality constraints
For improved efficiency, Mosek is solving the dual problem.
------------------------------------------------------------
MOSEK Version 8.0.0.60 (Build date: 2017-3-1 13:09:33)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 22
Cones : 18
Scalar variables : 105
Matrix variables : 0
Integer variables : 0
Optimizer started.
Conic interior-point optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 4
Eliminator terminated.
Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.02
Lin. dep. - number : 0
Presolve terminated. Time: 0.06
Interior-point optimizer terminated. Time: 0.11.
MOSEK DUAL INFEASIBILITY REPORT.
Problem status: The problem is dual infeasible
Optimizer terminated. Time: 0.30
Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -1.5773806652e-001 nrm: 1e+000 Viol. con: 0e+000 var: 0e+000 cones: 0e+000
Optimizer summary
Optimizer - time: 0.30
Interior-point - iterations : 0 time: 0.11
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
------------------------------------------------------------
Status: Infeasible
Optimal value (cvx_optval): -Inf
But when I replace the ‘-rel_entr(1, S_m_v_cat(:,k)+A_m_k_1(:,k))’ with the ‘log(S_m_v_cat(:,k)+A_m_k_1(:,k))’. The problem can be solved. It can be seen as below.
function [S_m_out, fval_2, flag_2] = S_m_calc_CVX(alpha, S_m_r, p_m, q_m_x, q_m_y)
%%参数定义
global P_charge;
global H;
global d_min;
global V_max;
global S;
global P_max;
global freq; %载波频率
global light_speed; %光速
global B;
global yeta_los;
global delta2;
global K;
global M;
global N;
global R_th;
global P_user;
%%变量初始化
A_m_k = zeros(N,M,K);
for m=1:1:M %%计算 c_m_k'(n)
for k=1:1:K
for n=1:1:N
dist2 = H*H + (q_m_x(n,m)-P_user(1,k))^2 + (q_m_y(n,m)-P_user(1,k))^2;
L_current = (4*pi*freq/light_speed)^2 * dist2 * yeta_los;
P_receive_current = p_m(n,m)/L_current;
P_receive_total = 0;
for m1=1:1:M;
dist2 = H*H + (q_m_x(n,m1)-P_user(1,k))^2 + (q_m_y(n,m1)-P_user(1,k))^2;
L_current = (4*pi*freq/light_speed)^2 * dist2 * yeta_los;
P_receive_temp = p_m(n,m1)/L_current;
P_receive_total = P_receive_total + P_receive_temp;
end
A_m_k(n,m,k) = P_receive_current/(P_receive_total - P_receive_current + B*delta2);
end
end
end
alpha_1 = reshape(alpha, N*M, K);
alpha_m_n_0 = sum(alpha_1,2);
A_m_k_1 = reshape(A_m_k, N*M, K);
S_m_r_0 = S_m_r;
for k=1:1:K-1
S_m_r_0 = cat(3,S_m_r_0,S_m_r);
end
alpha_m_n = alpha_m_n_0;
for k=1:1:K-1
alpha_m_n = cat(2,alpha_m_n,alpha_m_n_0);
end
S_m_r_1 = reshape(S_m_r_0, N*M, K);
% CC = -log(S_m_r_1)./log(2) + 1/log(2);
CC = log(S_m_r_1)./log(2);
cvx_begin
cvx_solver mosek
variable S_m_v(N,M)
expression S_m_v_cat
S_m_v_cat = S_m_v;
for k=1:1:K-1 %扩展S_m_r的维度到N*M*K,维度统一方便后续CVX计算
S_m_v_cat = cat(3,S_m_v_cat,S_m_v);
end
S_m_v_cat = reshape(S_m_v_cat, N*M, K);
maximize(sum(sum(alpha_1.*(log(S_m_v_cat+A_m_k_1)./log(2) - (CC + (S_m_v_cat - S_m_r_1) ./(S_m_r_1.*log(2))) ))).*B./(M*N));
subject to
for k=1:1:K
sum(alpha_1(:,k).*(log(S_m_v_cat(:,k)+A_m_k_1(:,k))./log(2) - (CC(:,k) + (S_m_v_cat(:,k) - S_m_r_1(:,k)) ./(S_m_r_1(:,k).*log(2))) )).*B./N >= R_th;
end
S_m_v >= sum(alpha,3);
S_m_v <= S;
cvx_end
S_m_out = reshape(S_m_v(:,1),N,M);
fval_2 = cvx_optval;
flag_2 = cvx_status;
where the output shows:
CVX Warning:
Models involving "log" or other functions in the log, exp, and entropy
family are solved using an experimental successive approximation method.
This method is slower and less reliable than the method CVX employs for
other models. Please see the section of the user's guide entitled
The successive approximation method
for more details about the approach, and for instructions on how to
suppress this warning message in the future.
Successive approximation method to be employed.
For improved efficiency, Mosek is solving the dual problem.
Mosek will be called several times to refine the solution.
Original size: 160 variables, 58 equality constraints
40 exponentials add 320 variables, 200 equality constraints
-----------------------------------------------------------------
Cones | Errors |
Mov/Act | Centering Exp cone Poly cone | Status
--------+---------------------------------+---------
30/ 30 | 3.215e-01 6.415e-03 0.000e+00 | Solved
20/ 20 | 2.404e-02 3.609e-05 0.000e+00 | Solved
8/ 8 | 3.421e-03 6.713e-07 0.000e+00 | Solved
0/ 0 | 4.105e-04 0.000e+00 0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +21449.8
I was confused by the problem. Can any one help me? Thanks in advance.