What is the correct form to express an SOCP problem in YALMIP?

[叶超天]

This is my problem, the variables are phi and psi.

phi=sdpvar(N,1); 

psi=sdpvar(N,1);

CS=[phi>0,psi>0];

CS=[CS,(GammaUCon+1)*abs(T_B1(1,:)*G_B*F_1(:,1))^2*(psi(1)^2)>=...

    GammaUCon*(norm(T_B1(1,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B1(1,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B1(1,:),2)^2)];

CS=[CS,(GammaUCon+1)*abs(T_B2(1,:)*G_B*F_2(:,1))^2*(psi(2)^2)>=...

    GammaUCon*(norm(T_B2(1,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B2(1,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B2(1,:),2)^2)];

CS=[CS,(GammaUCon+1)*abs(T_B2(2,:)*G_B*F_2(:,2))^2*(psi(3)^2)>=...

    GammaUCon*(norm(T_B2(2,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B2(2,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B2(2,:),2)^2)];

CS=[CS,(GammaUCon+1)*abs(T_B3(1,:)*G_B*F_3(:,1))^2*(psi(4)^2)>=...

    GammaUCon*(norm(T_B3(1,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B3(1,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B3(1,:),2)^2)];

CS=[CS,(GammaUCon+1)*abs(T_B3(2,:)*G_B*F_3(:,2))^2*(psi(5)^2)>=...

    GammaUCon*(norm(T_B3(2,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B3(2,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B3(2,:),2)^2)];

CS=[CS,(GammaUCon+1)*abs(T_B3(3,:)*G_B*F_3(:,3))^2*(psi(6)^2)>=...

    GammaUCon*(norm(T_B3(3,:)*G_B*F*diag(psi),'fro')^2+sigma_R^2*norm(T_B3(3,:)*G_B*F*inv(V),'fro')^2+sigma_B^2*norm(T_B3(3,:),2)^2)];

CS=[CS,norm(inv(H_B)*V*diag(phi),'fro')^2<=P_B];

CS=[CS,norm(inv(H_M)*V1*psi(1),'fro')^2<=P_M1];

CS=[CS,norm(inv(H_M)*V2*diag(psi(2:3)),'fro')^2<=P_M2];

CS=[CS,norm(inv(H_M)*V3*diag(psi(4:6)),'fro')^2<=P_M3];

CS=[CS,norm(F*diag(phi),'fro')^2+norm(F*diag(psi),'fro')^2+sigma_R^2*norm(F*inv(V),'fro')^2<=P_R];

Ipsi_1=diag([0;psi(2:6)]); % 

 Ipsi_2=diag([psi(1);0;0;psi(4:6)]);

 Ipsi_3=diag([psi(1:3);0;0;0]);

MT_1=[diag(phi) Ipsi_1];

MT_2=[diag(phi) Ipsi_2];

 MT_3=[diag(phi) Ipsi_3];

CS=[CS,(1+t_zero)*abs(T_M1(1,:)*G_M1*F_1(:,1))^2*(phi(1)^2)>=...

t_zero*(norm(T_M1(1,:)*G_M1*F*MT_1,'fro')^2+sigma_R^2*norm(T_M1(1,:)*G_M1*F*inv(V),'fro')^2+sigma_M^2*norm(T_M1(1,:),2)^2)];

CS=[CS,(1+t_zero)*abs(T_M2(1,:)*G_M2*F_2(:,1))^2*(phi(2)^2)>=...

t_zero*(norm(T_M2(1,:)*G_M2*F*MT_2,'fro')^2+sigma_R^2*norm(T_M2(1,:)*G_M2*F*inv(V),'fro')^2+sigma_M^2*norm(T_M2(1,:),2)^2)];

CS=[CS,(1+t_zero)*abs(T_M2(2,:)*G_M2*F_2(:,2))^2*(phi(3)^2)>=...

t_zero*(norm(T_M2(2,:)*G_M2*F*MT_2,'fro')^2+sigma_R^2*norm(T_M2(2,:)*G_M2*F*inv(V),'fro')^2+sigma_M^2*norm(T_M2(2,:),2)^2)];

CS=[CS,(1+t_zero)*abs(T_M3(1,:)*G_M3*F_3(:,1))^2*(phi(4)^2)>=...

t_zero*(norm(T_M3(1,:)*G_M3*F*MT_3,'fro')^2+sigma_R^2*norm(T_M3(1,:)*G_M3*F*inv(V),'fro')^2+sigma_M^2*norm(T_M3(1,:),2)^2)];

CS=[CS,(1+t_zero)*abs(T_M3(2,:)*G_M3*F_3(:,2))^2*(phi(5)^2)>=...

t_zero*(norm(T_M3(2,:)*G_M3*F*MT_3,'fro')^2+sigma_R^2*norm(T_M3(2,:)*G_M3*F*inv(V),'fro')^2+sigma_M^2*norm(T_M3(2,:),2)^2)];

CS=[CS,(1+t_zero)*abs(T_M3(3,:)*G_M3*F_3(:,3))^2*(phi(6)^2)>=...

t_zero*(norm(T_M3(3,:)*G_M3*F*MT_3,'fro')^2+sigma_R^2*norm(T_M3(3,:)*G_M3*F*inv(V),'fro')^2+sigma_M^2*norm(T_M3(3,:),2)^2)];

sol=solvesdp(CS,[],ops);

I know the above problem is a socp problem, but not in standard form. Will yalmip accept this form? I think the answer is no, cause the answer to the above problem is always

imfeasible, and the answers to phi and psi are both all 0 vectors. I have also tried to put a sqrt function in both sides, the answer is also imfeasible and always all 0 vectors.

The value I gave to the problem for verification ensures that the problem is feasible. Then a contradiction occurs. So could you tell me how to make yalmip accept this problem and give

the right answer.

Thank you very much.

[Johan Löfberg]

The fact that the problem has been solved, albeit proven to be infeasible, indicates that YALMIP accepts it. Infeasibility has nothing to do whether it is an SOCP or not. norm(x) <= -1 is an SOCP, but trivially infeasible.

If the problem is infeasible, the values you extract has no meaning as those values just happen to be some stuff returned by the solver when it gave up and claimed infeasibility.

Without reproducible code, no more can be said.

[叶超天]

Dear  Johan Löfberg:

       It is a great pleasure to hear from you. Now I have a  reproducible code. 

       

       First, I get phiC, psiC.  (line 66&67)

      

       The above phiC and psiC are a feasible  solution to the socp problem

       (phiC and psiC are diagonal matrix, their diagonal elements are used as phi and psi)

       all the constraints are satisfied when I substitute phiC and psiC in the problem except that

       line 181

       norm(inv(H_M)*V3*b3*eye(3),'fro')^2

       ans =

       300.0000

       P_M3

       P_M3 =

       300

       norm(inv(H_M)*V3*b3*eye(3),'fro')^2<=P_M3

       ans =0

       I think in fact this constraint is also satisfied.

       But I still can not use solvesdp to prove that the problem is feasible. 

       Please help me.

       Thanks a lot!

[叶超天]

       all the constraints are satisfied when I substitute phiC and psiC in the problem except that

       line 152

[Johan Löfberg]

First, there should be a lot of warnings when you run your code

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code

Second, have you noticed which solver was selected to solve your problem? It can not be written as a linear SOCP problem since you have quadratics on the wrong side of the inequalities in the norm constraints. Instead, YALMIP writes the problem as a general nonconvex quadratic program and, on my machine, calls the genreal nonlinear solver fmincon, which easily finds a feasible solution sometimes, and fails sometimes

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    0.000000e+00    9.966e+01    0.000e+00
    1       2    0.000000e+00    9.698e+01    0.000e+00    2.264e+00
    2       4    0.000000e+00    7.395e+01    0.000e+00    8.678e+00
    3       6    0.000000e+00    0.000e+00    0.000e+00    1.385e+01

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in
feasible directions, to within the selected value of the function tolerance,
and constraints are satisfied to within the selected value of the constraint tolerance.

<stopping criteria details>

sol =

    yalmiptime: 0.3900
    solvertime: 0.0310
          info: 'Successfully solved (FMINCON)'
       problem: 0
>>
>> check(CS)
 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|    ID|                           Constraint|   Primal residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|    #1|               Elementwise inequality|             2.975|
|    #2|               Elementwise inequality|            7.5579|
|    #3|   Elementwise inequality (quadratic)|           0.77358|
|    #4|   Elementwise inequality (quadratic)|           0.37723|
|    #5|   Elementwise inequality (quadratic)|           0.66346|
|    #6|   Elementwise inequality (quadratic)|           0.41113|
|    #7|   Elementwise inequality (quadratic)|           0.23162|
|    #8|   Elementwise inequality (quadratic)|            0.3089|
|    #9|   Elementwise inequality (quadratic)|           171.079|
|   #10|   Elementwise inequality (quadratic)|           85.1192|
|   #11|   Elementwise inequality (quadratic)|           14.8055|
|   #12|   Elementwise inequality (quadratic)|           183.163|
|   #13|   Elementwise inequality (quadratic)|          657.6619|
|   #14|   Elementwise inequality (quadratic)|           0.97976|
|   #15|   Elementwise inequality (quadratic)|           0.72895|
|   #16|   Elementwise inequality (quadratic)|           0.95692|
|   #17|   Elementwise inequality (quadratic)|           0.98143|
|   #18|   Elementwise inequality (quadratic)|           0.82421|
|   #19|   Elementwise inequality (quadratic)|           0.12178|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 

[Johan Löfberg]

You are comparing two numbers which are 300 when truncated in display. Up to the solvers precision, the constraint is satisfied. You cannot expect a tolerance any better than, say, 1e-8

>> A = 300

A =

   300

>> B = 300-1e-15

B =

   300

>> A<B

ans =

     0

[Johan Löfberg]

I see now that you are accidentally destroying an originally convex form.

The SOCP

norm([a;b]) <= c

should not be written as a^2 + b^2 <= c^2 which essentially is what you are doing

So, norm(a,'fro')^2 + norm(b,'fro')^2 <= c^2, I would write as

sdpvar s t
norm(a,'fro') <= t
norm(b,'fro') <= s
norm([t;s]) <= c % t^2+s^2 <= c^2

[叶超天]

That may be the problem. I will use the form you presented to try again.

In fact, at first, I use sqrt() at both sides, so at one side it is a socp form, but at the other side it's sqrt(norm()^2+norm^2).

It seems that yalmip can not distinguish this form as a standard socp either, and the answer is also imfeasible. 

Thank you very much.

[Johan Löfberg]

You absolutely don't want to do sqrt(norm()^2+norm^2) <= norm() either, as that obviously violates convexity propagation rules (norm is convex, hence norm() > something is not necessarily convex), and of course, it cannot be written as an SOCP directly from that form

[叶超天]

Then if I convert the problem into a standard socp form manually, ||Ax+b||<=cTx+d, the problem will sure be solved, is that right?

[Johan Löfberg]

Of course, but you don't have to go that far. As long as you don't square the norms on the wrong side of the constraints etc, there is no problem

norm(x,2)^2 <= 1 is ok, but is a special case, as norm(x)^2 immediately is replaced with x'*x, hence it simply says x'*x <= 1

norm(x,1)^2 <= 1 is not good, as squaring a general convex operator is not necessarily something which preserves convexity (it is convex though, but YALMIP doesn't look for this special case)

norm(x,2)^2 <= y^2 is not good, as you have a convex operator (y^2) on the wrong side of the inequality

norm(x,2)^2 <=norm(y) is not good, as you have a convex operator (norm(y)) on the wrong side of the inequality

Just don't apply any operators on the norm expressions, and nothing can go wrong, as long as you understand on which side you can have them

and as I said before

norm(x,1)^2 + norm(x,2)^2 <= z^2

is

s^2 + t^2 <= z^2, norm(x,1)<=s, norm(x,2)<= t

which is

norm([s;t]) <= z (assuming z>=0), norm(x,1)<=s, norm(x,2)<= t

which satisfies all convexity propagation rules

[叶超天]

At last, the program works. I'm just going to convert it to the standard form manually. Your method save me a lot of time.  I feel very happy now. Many thanks.

I suggest  you  consider add a new general standard SOCP form which is  easy to use in yalmip, such that we do not have to add new variables. The program can convert it automatically. Just a suggestion.

Since you are professional in optimization, I think it is a good chance to ask you a few more questions and learn from you.

1 I 'm studying communication problems and often  encounter problems which involve complex matrices. I wonder if there is a way to decide whether an objective function 

is convex or not. Is there some hidden rules, experience, or techniques. Sometimes, somebody can tell whether the problem is convex or not immediately, what's the secret? Are  there any  books or materials that you recommend?

As far as I know, the following  methods can be used to decide convexity:

1) the hessian matrix 

2) the propagation rules you mentioned (but I'm not very good at using this method)

3) two unknown matrix variables multiplied together is not convex.

4) if the problem is not convex in some of the variables, then it is not convex in all the variables. (is this right?)

5) For a problem which has vector variables, if the low dimensional version of the problem is nonconvex, then the original problem is nonconvex.

What is your comment?

2 I have seen some similar "substitution technique" in Stephen Boyd's book, by which I mean 

sdpvar s t
norm(a,'fro') <= t
norm(b,'fro') <= s
norm([t;s]) <= c % t^2+s^2 <= c^2 

I did realize to use this method. This is really important.

 Could you tell me when to use that and more examples, or related materials.

[叶超天]

I did realize to use this method. -> I did NOT realize to use this method. 

[Johan Löfberg]

YALMIP already does a lot of reformulations and introduces a lot of variables behind your back, but you have one which cannot be dealt with.

norm(x)^2 <= y^2

This cannot be reformulated automatically as it is a nonconvex set ((1,1) is feasible as is (1,-1), but not (1,0) inbetween). The reason I could give you a work-around is that I guessed that you had the constraint y>=0 also, although not explicitly stated.

Regarding your questions, YALMIP convexity analysis is entirely based on simple propagation rules (along with the basic notation that the set f(x)<=g(x) is convex if f(x) is convex and g(x) is concave etc)
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.GraphRepresentations

Beyond that, once you've read, e.g., Boyd's book, it is all about practice and experience.