CVX using MOSEK bug


#1
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.


(Mark L. Stone) #2

You haven’t provided reproducible problems.

Try removing the objective e function for both formulations and see what happens. As was suggested to you at Failed using CVX. Please help , perhaps the scaling is poor. But the chances might be better to get a correct assessment of feasibility if there is not a possibly poorly scaled objective function to complicate things.


#3

Dear Mark
I will try your way and carefully check the scalling of my model. I hope this will work. Again, many thanks for your help.