Unable to solve Convex programming problem

I have the following data and code for the convex problem.

mean_c = [5; 15; 18 ; 26];
mean_d1 = [25 ; 21 ; 8 ; 16];
mean_d2 = [20 ; 25 ; 31 ; 33];
covariance_c_sqrt = [1.6424 0.6645 0.5900 0.8558;
0.6645 1.4634 0.9441 0.6344;
0.5900 0.9441 1.4894 0.5574;
0.8558 0.6344 0.5574 1.2512];
covariance_d1_sqrt = [1.5406 1.1921 1.2443 0.7342;
1.1921 1.6658 0.7972 0.7327;
1.2443 0.7972 1.7399 1.2046;
0.7342 0.7327 1.2046 1.2911];
covariance_d2_sqrt = [1.9069 1.3501 1.0056 0.2173;
1.3501 1.4972 0.6666 0.7804;
1.0056 0.6666 0.7978 0.4356;
0.2173 0.7804 0.4356 1.8474];
p0 = 0.9;
p1 = 0.9;
a = 0.7;
P1 = [0.5 0.5 ; 0.3 0.7];
P2 = [0.6 0.4 ; 0.1 0.9];
A1 = eye(2) - a.*P1;
A2 = eye(2) - a.*P2;
A = [A1 ; A2];
b = [0.5 0.5];
ub1 = 30;
ub2 = 40;
cvx_begin
variables rho(4,1) y(2,1) z(1)
minimize (z)
rho’*mean_c + sqrtm( p0/(1-p0) )norm(covariance_c_sqrtrho) <= z;
rho’*mean_d1 + sqrtm( (p1^y(1))/(1-p1^y(1)) )*V1 <= ub1;
rho’*mean_d2 + sqrtm( (p1^y(2))/(1-p1^y(2)) )*V2 <= ub2;
y(1) + y(2) == 1;
y >= 0;
rho >= 0;
rho’*A == (1-a)*b;
cvx_end

I get the output

Check for incorrect argument data type or missing argument in call to function ‘vec’.

Error in cvx/sparsify>replcols (line 142)
cvx___.readonly( ndxs ) = vec( cvx_readlevel( bN ) );

Error in cvx/sparsify (line 105)
[ x, forms, repls ] = replcols( x, tt, ‘none’, forms, repls, isobj );

Error in cvx/exp (line 68)
xt = sparsify( xt, ‘exponential’ );

Error in .^ (line 121)
cvx_optval = exp( log( cvx_constant( xt ) ) .* yt );

Error in ^ (line 9)
z = power( x, y );

Could you please suggest how do I transform it, I used log/square on both sides, still same output.
Also should I change it to cvxquad instead of cvx?

From DCP ruleset, I wrote the problem again as
cvx_begin
cvx_solver gurobi
variables rho(4,1) y(2,1) z
minimize (z)

          rho'*mean_c + sqrt( p0/(1-p0) )*norm( covariance_c_sqrt*rho ) <= z,    
         rho'*mean_d1 + sqrt( inv_pos( inv_pos( pow_p(p1,y(1)) ) - 1 ) )*V1 <= ub1;
          rho'*mean_d2 + sqrt( inv_pos( inv_pos( pow_p(p1,y(2)) ) - 1 ) )*V2 <= ub2;
          y(1) + y(2) == 1,
          y >= 0;
          rho >= 0;
          rho'*A == (1-a)*b;

cvx_end

But this gives infeasibility/unboundedness irrespective of the values I fix.
Maybe it is due to the 2nd and 3rd constraints. I am unable to simplify it any further.

pow_p(p1,y(1)) should not be allowed by CVX, because the 2nd argument is supposed to be a constant. However, CVX allows it and considers it to be a constant. And therefore it allows the iterated inv_pos which would be disallowed, but for CVX’s improperly thinking pow_p(p1,y(1)) is a constant. Not to mention the sqrt would then also be disallowed.

So whatever CVX is sending to the solver based on the 2nd and 3rd constraints is meaningless garbage. Therefore, any results produced by the solver and CVX from your program are meaningless garbage.

You should start with a CONVEX optimization problem. Then input it to CVX according to its DCP rules. Your current program violates CVX rules like it’s going out of style, and you’re rushing to cram as many rules violations into your constraints as you can.

help pow_p
pow_p Positive branch of the power function.
pow_p(X,P) computes a convex or concave branch of the power function:
P < 0: pow_p(X,P) = X.^P if X > 0, +Inf otherwise
0 <= P < 1: pow_p(X,P) = X.^P if X >= 0, -Inf otherwise
1 <= P : pow_p(X,P) = X.^P if X >= 0, +Inf otherwise
Both P and X must be real.

Disciplined convex programming information:
    The geometry of pow_p(X,P) depends on the precise value of P,
    which must be a real constant:
             P < 0: convex  and nonincreasing; X must be concave.
        0 <= P < 1: concave and nondecreasing; X must be concave.
        1 <= P    : convex  and nonmonotonic;  X must be affine.
    In all cases, X must be real.

Other functions named pow_p

help inv_pos
inv_pos Reciprocal of a positive quantity.
inv_pos(X) returns 1./X if X is positive, and +Inf otherwise.
X must be real.

 For matrices and N-D arrays, the function is applied to each element.

  Disciplined convex programming information:
      inv_pos is convex and nonincreasing; therefore, when used in CVX
      specifications, its argument must be concave (or affine).

Other functions named inv_pos

help cvx/sqrt
Discipined convex programming information for sqrt:
sqrt(X) is log-concave and nondecreasing in X. Therefore, when used
in DCPs, X must be concave (or affine).

Disciplined geometric programming information for sqrt:
   sqrt(X) is log-log-affine and nondecreasing in X. Therefore, when
   used in DGPs, X may be log-affine, log-convex, or log-concave.

Thank you for your reply.
I rewrote the code and used two solvers. The MOSEK solver gives the value while SDPT3 does not.

cvx_begin
cvx_solver mosek
variables rho(4,1) y(2,1) x(2,1) z(1)
minimize (z)

          rho'*mean_c + sqrt( p0/(1-p0) )*norm( covariance_c_sqrt*rho ) <= z, 

          0.5*y(1)*log(p1) - 0.5*log(1-x(1)) + log(V1) <= log(ub1 - rho'*mean_d1);
          0.5*y(2)*log(p1) - 0.5*log(1-x(2)) + log(V2) <= log(ub2 - rho'*mean_d2);

          y(1)*log(p1) <= log(x(1));
          y(2)*log(p1) <= log(x(2));
           
          
          y(1) + y(2) == 1,
          y >= 0;
          rho >= 0;
          rho'*A == (1-a)*b;

cvx_end

Calling Mosek 9.1.9: 31 variables, 10 equality constraints
For improved efficiency, Mosek is solving the dual problem.

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 10
Cones : 7
Scalar variables : 31
Matrix variables : 0
Integer variables : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 2
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 10
Cones : 7
Scalar variables : 31
Matrix variables : 0
Integer variables : 0

Optimizer - threads : 12
Optimizer - solved problem : the primal
Optimizer - Constraints : 7
Optimizer - Cones : 7
Optimizer - Scalar variables : 29 conic : 23
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 18 after factor : 19
Factor - dense dim. : 0 flops : 2.87e+02
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 6.9e+00 5.0e+00 1.9e+01 0.00e+00 1.481043233e+01 -3.000000000e+00 1.0e+00 0.00
1 2.6e+00 1.9e+00 6.7e+00 -2.33e-01 2.637733647e+00 -7.517742550e+00 3.9e-01 0.00
2 4.6e-01 3.3e-01 6.1e-01 4.16e-01 -1.204365302e+01 -1.424654918e+01 6.7e-02 0.00
3 9.1e-02 6.6e-02 5.1e-02 9.39e-01 -1.559793015e+01 -1.604148726e+01 1.3e-02 0.00
4 3.4e-02 2.5e-02 1.5e-02 8.26e-01 -1.611210768e+01 -1.630711868e+01 5.0e-03 0.00
5 7.5e-03 5.4e-03 2.1e-03 7.11e-01 -1.640435272e+01 -1.645459867e+01 1.1e-03 0.00
6 1.6e-03 1.1e-03 6.8e-04 1.85e-02 -1.649961733e+01 -1.651707211e+01 2.3e-04 0.00
7 7.1e-04 5.1e-04 3.2e-04 -4.98e-01 -1.669840670e+01 -1.670145038e+01 1.0e-04 0.00
8 2.6e-04 1.9e-04 8.3e-05 4.80e-01 -1.695143068e+01 -1.695115939e+01 3.7e-05 0.00
9 4.3e-05 3.1e-05 6.0e-06 7.66e-01 -1.711252919e+01 -1.711249572e+01 6.3e-06 0.00
10 1.4e-05 1.0e-05 1.2e-06 8.93e-01 -1.714181523e+01 -1.714173518e+01 2.0e-06 0.00
11 2.1e-07 1.6e-07 2.2e-09 9.88e-01 -1.715476387e+01 -1.715476358e+01 3.1e-08 0.00
12 5.3e-09 3.9e-09 8.5e-12 1.00e+00 -1.715497572e+01 -1.715497572e+01 7.8e-10 0.00
Optimizer terminated. Time: 0.00

Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: -1.7154975718e+01 nrm: 1e+02 Viol. con: 5e-07 var: 0e+00 cones: 0e+00
Dual. obj: -1.7154975907e+01 nrm: 6e+00 Viol. con: 0e+00 var: 9e-08 cones: 0e+00
Optimizer summary
Optimizer - time: 0.00
Interior-point - iterations : 12 time: 0.00
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: Solved
Optimal value (cvx_optval): +17.155

cvx_begin
cvx_solver sdpt3
variables rho(4,1) y(2,1) x(2,1) z(1)
minimize (z)

          rho'*mean_c + sqrt( p0/(1-p0) )*norm( covariance_c_sqrt*rho ) <= z, 

          0.5*y(1)*log(p1) - 0.5*log(1-x(1)) + log(V1) <= log(ub1 - rho'*mean_d1);
          0.5*y(2)*log(p1) - 0.5*log(1-x(2)) + log(V2) <= log(ub2 - rho'*mean_d2);

          y(1)*log(p1) <= log(x(1));
          y(2)*log(p1) <= log(x(2));
           
          
          y(1) + y(2) == 1,
          y >= 0;
          rho >= 0;
          rho'*A == (1-a)*b;

cvx_end

Successive approximation method to be employed.
For improved efficiency, SDPT3 is solving the dual problem.
SDPT3 will be called several times to refine the solution.
Original size: 31 variables, 10 equality constraints
6 exponentials add 48 variables, 30 equality constraints

Cones | Errors |
Mov/Act | Centering Exp cone Poly cone | Status
--------±--------------------------------±--------
Check for incorrect argument data type or missing argument in call to function ‘vec’.

Error in cvxprob/solve (line 250)
Anew2 = Anew * diag(sparse(vec(amult(orow,:))));

Error in cvx_end (line 88)
solve( prob );

As suggested in here, I have CVX directories above that of YALMIP.
Could you please suggest if there is any problem with my code or else how do I remove the error?

Mosek solved your problem, so maybe you should be content with that. It is more reliable than using CVX’s Successive approximation method. However, the Successive approximation method never even had a chance to run.

YALMIP and CVX have different vec. Apparently, that is why you had a problem. I suggest having all CVX directories ahead of YALMIP in the MATLAB path. Perhaps that will avoid the conflict you experienced. If you do that, you should remove the "useless’ CVX functions, such as logdet, which do nothing other than tell the user to use the version with _, e.g., log_det, but which have the same name as the YALMIP version, and hence can cause YALMIP to fail, because CVX"s useless logdet would be used instead of YALMIP’s logdet. So by eliminating the useless CVX functions, you use log_det with CVX, and use logdet with YALMIP. Similarly with geomean.

You can use
which -all vec
to see all the versions of vec in the MATLAB path, and their order in the path.

which -all logdet
etc.

1 Like

Noted. Thank you so much!