Search code examples
pythonoptimizationcvxpyconvex-optimization

CVXPY: DCPError: Problem does not follow DCP rules


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.sstatic.net/IQpxh.png

Solution

  • 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)))