2

I am trying to code a problem solves B(2,1) under LMI constraints.

R(2,1)=R0(2,1)+H(2,2)*B(2,1)

Vc is a scalar variable

It keeps getting

> "DCPError: Problem does not follow DCP rules."

    import numpy as np
    import cvxpy as cp

    H = np.random.rand(2,2)
    R0 = np.random.rand(2,1)

    B=cp.Variable((2,1), complex=True)  
    Rf=cp.diag(R0+H*B)




    RRf=cp.real(Rf)
    IRf=cp.imag(Rf)


    Vc=cp.Variable()
    Vc2= (Vc**2)


    z=np.zeros((Rf.shape[0],Rf.shape[1]))
    I=np.eye(Rf.shape[0])
    objective3=cp.Minimize(Vc2)

    LMI =cp.bmat( [   
                            [Vc2*I,        RRf,    z,        -IRf],
                            [RRf,          I,      IRf,          z],
                            [z,            IRf,    Vc2*I,        RRf],
                            [-IRf,         z,       RRf,          I]      
                                                                                ]) 
    const1 = LMI  >=0
    const2 = Vc   >=0        

    prob=cp.Problem(objective3,[const1,const2])
    print(prob.is_dcp())   




  [1]: https://i.stack.imgur.com/IQpxh.png
Maged Eltorkoman
  • 67
  • 1
  • 7
  • Provided that the objective function is DCP =True and ----- LMI block is DCP=True ----- const1 DCP is = False though ?! – Maged Eltorkoman Apr 27 '20 at 07:16
  • You have the square of a variable appearing inside the matrix. That is not DCP compliant. – Michal Adamaszek Apr 28 '20 at 16:38
  • Thanks for replying! @MichalAdamaszek, I replaced the variable Vc instead of the its square in the LMI block and it returned a valid DCP problem .. unfortunately, the following error appeared: **Anaconda3\lib\site-packages\cvxpy\reductions\complex2real\complex2real.py", line 146, in canonicalize_expr assert all(v is None for v in imag_args) AssertionError** – Maged Eltorkoman Apr 28 '20 at 18:37
  • As far as I know to get an LMI you should write ``LMI>>0`` and not ``LMI>=0``. What you did now is just say that all entries are nonnegative, which may explain the error (imaginary parts are nonzero). – Michal Adamaszek Apr 28 '20 at 19:17
  • @MichalAdamaszek I tried `LMI>>0` but it gives the same **AssertionError**. Can it be because I am using a complex variable?
    I am really stuck with this issue for weeks by now. I am a mechanical engineer not a programmer but this code would be a great help at work.     – Maged Eltorkoman Apr 28 '20 at 20:17
  • I ran your code now and you are right that something is wrong. Perhaps it is a bug in the real/imag operators for complex variables that should be reported to cvxpy. Anyway, it worked for me when I defined the real and imaginary parts separately: ``BR=cp.Variable((2,1)) BI=cp.Variable((2,1)) RRf=cp.diag(R0+H@BR) IRf=cp.diag(H@BI)`` – Michal Adamaszek Apr 28 '20 at 20:49
  • @MichalAdamaszek You are absolutely right the code works that way but the math doesn't!! But your suggestion was a great help the right way ... THANKS!!! – Maged Eltorkoman Apr 29 '20 at 16:00
  • You were right in the first place . It should be `LMI>>0` not `LMI >=0` but when you put it `>>` the standard SCS solver doesn't handle the matrix size. So I downloaded the CVXOPT solver and it worked fine till now. I hope no more bugs so I can complete the whole thing. – Maged Eltorkoman Apr 29 '20 at 16:05

1 Answers1

1

With @MichalAdamaszek 's help the next code works.
The problem was that CVXPY is not able to handle .real and .imag function inside the constraints.
So the manipulation needed was to break down the complex varibale B into two real variables then combine them after the .solve usingB=BR.value+1j*BI.value
The other mistake in the question was putting the constraint as LMI>=0. For SDP LMI>>0 should be used. The last thing was to use CVXOPT solver instead the standardSCS as it can't handle more than 2x2 matrices. The code proves to be mathmatically correct as it always minimize the residual function

R(2,1)=R0(2,1)+H(2,2)*B(2,1)

print('The residule',abs(R0+np.matmul(H,B))) approaches 0 every run time.
The correct code: 

import numpy as np
import cvxpy as cp

H = np.random.rand(2,2)
R0 = np.random.rand(2,1)

BR=cp.Variable((2,1))
BI=cp.Variable((2,1))  




RRf=cp.diag((np.real(R0)+np.real(H)@BR-np.imag(H)@BI))
IRf=cp.diag((np.imag(R0)+np.imag(H)@BR+np.real(H)@BI))

Vc2=cp.Variable()


z=np.zeros((RRf.shape[0],RRf.shape[1]))
I=np.eye(RRf.shape[0])
objective3=cp.Minimize(Vc2)


LMI =cp.bmat( [   
                        [Vc2*I,        RRf,    z,        -IRf],
                        [RRf,          I,      IRf,          z],
                        [z,            IRf,    Vc2*I,        RRf],
                        [-IRf,         z,       RRf,          I]      
                                                                            ]) 
const1 = LMI  >>0



prob=cp.Problem(objective3,[const1])
prob.solve(solver=cp.CVXOPT, kktsolver=cp.ROBUST_KKTSOLVER)
B=BR.value+1j*BI.value

print(abs(B),Vc2.value)
print('The residule',abs(R0+np.matmul(H,B)))
Maged Eltorkoman
  • 67
  • 1
  • 7