Why CVX completely ignores a constraint or violates it?

Mobin · April 30, 2020, 7:52am

Hello everybody

I am trying to solve a energy management problem using CVX. There are two categories of equality conditions including active power balance and reactive power balance. In addition, there are the inequality constraints, which constrain the power of generators and loads within acceptable ranges. I use CVX 2.1 and CVX 3 beta.

Using CVX 2.1: If I consider all equality and inequality constraints, the equality power balance constraints are not satisfied. When I remove the inequality constrains, the equality constraints is able to be satisfied. I don’t know why CVX ignores the equality constraints in the first situation.
Using CVX 3 beta: the following error will be appeared with mentioning to the cvx_end line in the main program:
>>> Subscript indices must either be real positive integers or logicals.<<<

The code is as follows. In order to see the violation, please check trace(C_p(:,:,n)*M(:,:,t))==P_inj(n,t) and trace(C_q(:,:,n)*M(:,:,t))==Q_inj(n,t) for example for n=4,5, or 6 and some times.

I will be so grateful if you guide me in this case.

clc
clear
close all

TM = 1*24;    
t0 = 1:TM;
P_d = 1*[50 45 43 42 43 50 64 65 59 59 58 59 59 60 59 61 63 65 93 100 96 92 74 50];
Q_d = P_d*tan(acos(0.8));
P_d1 = 0.4*P_d; 
P_d2 = 0.3*P_d; 
P_d3 = 0.3*P_d;
Q_d1 = 0.4*Q_d;
Q_d2 = 0.3*Q_d;
Q_d3 = 0.3*Q_d; 
P_cr = [12.5*ones(1,16) 25*ones(1,8)];
Q_cr = P_cr*tan(acos(0.8));
E_i = [250 500 750]; 
f_c = 0.01*[3 3 2.8 2.6 2.5 2.6 2.8 2.9 3.3 3.35 3.4 3.5 3.55 3.5 3.4 3.3 3.2 3.2 3.5 4 4.45 3.7 3.5 3.4];
f_r = 0.9*f_c; 
GM = [1 1 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 1 1 1]; 
IM=(GM==0); 
NoV = 60; 
n_L = 15; 
S_Tmax = 500; 
S_gmax= 50;  
P_gmax = 0.8*S_gmax; 
Q_gmax = sin(acos(0.8))*S_gmax; 
G = 3; 
CL = 3; 
OL = 3; 
CR = 1;
fi = 0.0382; 
fd = 0.0427; 
fcr = 0.0441;
Ps = 0.3*[0 0 0 0 0 0 15 37 70 100 120 140 150 150 148 130 113 78 56 30 0 0 0 0];
Pw = 0.3*[150 150 150 130 128 127 120 120 122 125 125 126 125 125 123 80 63 60 55 60 75 82 105 130]; 
f_g = [0.0158 0.0135 0.0113];
h_g = [.0079 .0067 .0056]; 
lambda_g = [2 2]; 
r_g = (S_gmax*0.9/(0.01*50))*ones(1,G); 
Df_min = 0.5; %Hz
Df_max = -0.5; %Hz

Sum_TraceP = 0;
Sum_TraceQ = 0;

% ===================== Cost Function =======================

cvx_begin
% cvx_solver sedumi
variable P(NoV,length(t0))
expressions P_inj(18,24) Q_inj(18,24) P_hat_grid(1,24) P_ord1(1,24) P_ord2(1,24) P_ord3(1,24) P_critic(1,24) P_pv1(1,24) P_pv2(1,24) P_pv3(1,24) P_wind1(1,24) P_wind2(1,24) P_wind3(1,24) M(31,31,24) I_Ldot(31,1)

F1 = 0;
F2 = 0;
F3 = 0;
F4 = 0;
F5 = 0;
F6 = 0;

P_hat_grid = max(P(4,:)+P(53,:)+P(54,:),zeros(1,TM));
for t=1:TM

        F1 = F1 + (f_c(t)-f_r(t))*P_hat_grid(t);
end

for t=1:TM
        F2 = F2 + (f_r(t))*(P(4,t)+P(53,t)+P(54,t));
end

for i=1:G
    for t=1:TM
        F3 = F3 + f_g(i)*P(i,t) + h_g(i)*(P(i,t)^2);
    end
end

for i=5:CL+4
    for t=1:TM
       F4=F4 + fi*P(i,t); 
    end
end

for d=8:OL+7
     for t=1:TM
       F5 = F5+ fd*P(d,t);
     end
end

  for t=1:TM;
   F6 = F6+ fcr*P(11,t);
  end
 
F_total = F1 + F2 + F3 + F4 + F5 + F6; % Total Cost
minimize F_total
subject to

% ============ Operational Constraints =================
P(12,1:3)==[0 0 0];
P(12,20:24)==[0 0 0 0 0];
P(13,1:4)==[0 0 0 0];
P(13,21:24)==[0 0 0 0];
P(14,1:5)==[0 0 0 0 0];
P(14,22:24)==[0 0 0];

sum(P(12,:))== E_i(1) - sum(P(5,:));
sum(P(13,:))== E_i(2) - sum(P(6,:));
sum(P(14,:))== E_i(3) - sum(P(7,:));

P(15,1) == (100/100)* (E_i(1)); 
SE_min = 0.05*(E_i(1)); 
SE_max = (200/100)*(E_i(1)); 
Kesi = 0.9;

for t=1:TM-1
P(15,t+1) == Kesi*P(15,t)+P(12,t+1);
end

for t=1:TM
    SE_min<=P(15,t)<=SE_max;
end

Pi_max = 50; 
Pi_min1 = -SE_max; 

for t=1:TM
    Pi_min1<=P(12,t)<=Pi_max;
    0<=P(13,t)<=Pi_max;      
    0<=P(14,t)<=Pi_max;      
end

for i=5:CL+4
    0<=sum(P(i,:))<=E_i(i-4);
end

for t=1:TM
      0<=P(8,t)<=P_d1(t);
      0<=P(9,t)<=P_d2(t);
      0<=P(10,t)<=P_d3(t);
      0<=P(11,t)<=P_cr(t);
end

P_ord1 = P_d1-P(8,:);
P_ord2 = P_d2-P(9,:); 
P_ord3 = P_d3-P(10,:);
P_critic = P_cr-P(11,:);

P_pv1 = Ps/3-P(16,:); 
P_pv2 = Ps/3-P(57,:);
P_pv3 = Ps/3-P(58,:);

for t=1:TM
    0<=P(16,t)<=Ps(t)/3;
    0<=P(57,t)<=Ps(t)/3;
    0<=P(58,t)<=Ps(t)/3;
end

P_wind1 = Pw/3-P(17,:); % Produced active power of the wind generator 
P_wind2 = Pw/3-P(59,:);
P_wind3 = Pw/3-P(60,:);

for t=1:TM
    0<=P(17,t)<= Pw(t)/3;
    0<=P(59,t)<= Pw(t)/3;
    0<=P(60,t)<= Pw(t)/3;
end

zeta = 0.17; % Ref [10]
for t=1:TM
    (P(4,t)^2)+(P(18,t)^2)<= (1/9)*S_Tmax^2;
    (P(53,t)^2)+(P(55,t)^2)<= (1/9)*S_Tmax^2;
    (P(54,t)^2)+(P(56,t)^2)<= (1/9)*S_Tmax^2;
end

for i=19:G+18
for t=1:TM
    0<=P(i-18,t)<=P_gmax;
    0<=P(i,t)<=Q_gmax;
    (P(i,t)^2)+(P(i-18,t)^2)<= S_gmax^2;
end
end

for t=1:TM
    GM(t)*(P(16,t))==0;
    GM(t)*(P(57,t))==0;
    GM(t)*(P(58,t))==0;
    GM(t)*(P(17,t))==0;
    GM(t)*(P(59,t))==0;
    GM(t)*(P(60,t))==0;
    GM(t)*P(22,t)==0;
end

for t=1:TM           
    IM(t)*P(22,t)==0; % If a secondary frequency control is used
end

for t=1:TM            
    IM(t)*(P(4,t))==0;     
    IM(t)*(P(53,t))==0;
    IM(t)*(P(54,t))==0;     
    IM(t)*(P(18,t))==0;
    IM(t)*(P(55,t))==0;     
    IM(t)*(P(56,t))==0;
end

% Injected active and reactive powers to the buses
P_inj = [P(4,:);P(53,:);P(54,:);P(1,:);P(2,:);P(3,:);-P(12,:);-P(13,:);-P(14,:);P_pv1;P_pv2;P_pv3;-P_ord1;-P_ord2;-P_ord3-P_critic;P_wind1;P_wind2;P_wind3];
Q_inj = [P(18,:);P(55,:);P(56,:);P(19,:);P(20,:);P(21,:);zeros(1,24);zeros(1,24);zeros(1,24);zeros(1,24);zeros(1,24);zeros(1,24);-Q_d1;-Q_d2;-Q_d3-Q_cr;zeros(1,24);zeros(1,24);zeros(1,24)];

% Admitance matrix for the network
Z_se = (0.08+1j*0.05)*eye(3);
Y11 = inv(Z_se);
Y66 = inv(Z_se);
Y22 = 2*eye(3)/(Z_se);
Y33 = 2*eye(3)/(Z_se);
Y44 = 2*eye(3)/(Z_se);
Y55 = 2*eye(3)/(Z_se);
Y12 = -inv(Z_se);
Y21 = Y12;
Y23 = -inv(Z_se);
Y32 = Y23;
Y34 = -inv(Z_se);
Y43 = Y34;
Y45 = -inv(Z_se);
Y54 = Y45;
Y56 = -inv(Z_se);
Y65 = Y56;

Y = [Y11      Y12      zeros(3) zeros(3) zeros(3) zeros(3);...
     Y21      Y22      Y23      zeros(3) zeros(3) zeros(3);...
     zeros(3) Y32      Y33      Y34      zeros(3) zeros(3);...
     zeros(3) zeros(3) Y43      Y44      Y45      zeros(3);...
     zeros(3) zeros(3) zeros(3) Y54      Y55      Y56;...
     zeros(3) zeros(3) zeros(3) zeros(3) Y65      Y66]; 

Y1_11 = Y(1:3,1:3);
Y1_12 = Y(1:3,4:18);
Y1_21 = Y(4:18,1:3);
Y1_22 = Y(4:18,4:18);

A1_prim = [Y1_12/Y1_22;eye(n_L)];
A1_sec = [Y1_11-Y1_12*(Y1_22\Y1_21);zeros(n_L,3)];

B1_prim = [zeros(3,n_L);inv(Y1_22)];
B1_sec = [eye(3);-Y1_22\Y1_21];

V_grid = (380/sqrt(3))*[1;exp(2*pi*1j/3);exp(4*pi*1j/3)]; % Three-phase Grid voltages
V_nom = 380/sqrt(3);
V_min = 0.95*V_nom;
V_max = 1.05*V_nom;
THD_max = 0.05;
aT = zeros(18,15);
bT = zeros(18,15);
alpha = zeros(18,1);
betha = zeros(18,1);
ro = zeros(18,1);
Khe = zeros(18,1);
C_p = zeros(31,31,18);
C_q = zeros(31,31,18);
C_v1 = zeros(31,31,18);

load('PL0_WONetCon.mat')
load('QL0_WONetCon.mat')

for t=1:TM

I_Ldot = [P(23:37,t);P(38:52,t);1];

% Equilibrium point
S_L0 = PL0(:,t)+1j*QL0(:,t);
V_B0 = V_nom*[1;exp(2*pi*1j/3);exp(4*pi*1j/3);1;exp(2*pi*1j/3);exp(4*pi*1j/3);1;exp(2*pi*1j/3);exp(4*pi*1j/3);1;exp(2*pi*1j/3);exp(4*pi*1j/3);1;exp(2*pi*1j/3);exp(4*pi*1j/3)];
I_L0 = conj(S_L0)./conj(V_B0);
I_L0dot = [real(I_L0);imag(I_L0);0];
M(:,:,t) = I_L0dot*transpose(I_Ldot) + I_Ldot*transpose(I_L0dot) + I_L0dot*transpose(I_L0dot); % M = DM + M0
I1_dot = [real(I_L0);imag(I_L0)];
   
M(:,:,t)==semidefinite(2*n_L+1);
M(:,:,t)>=0;

for n=1:18
    
    aT(n,:) = A1_prim(n,:);
    bT(n,:) = B1_prim(n,:);
    alpha(n) = A1_sec(n,:)*V_grid;
    betha(n) = B1_sec(n,:)*V_grid;
    
    % Matrices for active powers of buses
    K_p = 0.5*((aT(n,:))'*bT(n,:)+(bT(n,:))'*aT(n,:));
    k_p = 0.5*(betha(n)*(aT(n,:))'+alpha(n)*(bT(n,:))');
    pay_p = real(conj(alpha(n))*betha(n));
    K_p2 = [real(K_p) -imag(K_p);imag(K_p) real(K_p)];
    k_p1 = [real(k_p);imag((k_p))]; % k_pdot
    C_p(:,:,n) = [K_p2 k_p1;transpose(k_p1) pay_p];
    
    
    % Matrices for reactive powers of buses
    Khe(n) = imag(conj(alpha(n))*betha(n));   
    q = (1/(2*1j))*(betha(n)*(aT(n,:))' - alpha(n)*(bT(n,:))');
    q_dot = [real(q);imag(q)];
    ro_q = Khe(n)^2 - (zeta^2)*pay_p^2;
    r_q = 2*(Khe(n)*q_dot - (zeta^2)*pay_p*k_p1) ;
    R_q = 4*(q_dot*transpose(q_dot) - (zeta^2)*k_p1*transpose(k_p1));
    C_q(:,:,n) = [R_q r_q;transpose(r_q) ro_q];
    
   
    % Matrices for voltage of buses
    S_V = (bT(n,:))'*bT(n,:);
    S_V2 = [real(S_V) -imag(S_V);imag(S_V) real(S_V)];
    s_v = betha(n)*(bT(n,:))';
    s_v1 = [real(s_v);imag(s_v)];
    sigma_v = (abs(betha(n)))^2;
    nu = 0.05; % Voltage tolerance
    C_v1(:,:,n) = [S_V2 s_v1;transpose(s_v1) sigma_v-((1+nu)^2)*V_nom^2];
    
    
    % >>>> Active power balance and reactive power balance constraints <<<
    trace(C_p(:,:,n)*M(:,:,t))==P_inj(n,t);
    trace(C_q(:,:,n)*M(:,:,t))==Q_inj(n,t);
    
    %  and voltage constraints
    trace(C_v1(:,:,n)*M(:,:,t))<=V_max^2;
    trace(C_v1(:,:,n)*M(:,:,t))>=V_min^2;
end

end

cvx_end

Michal_Adamaszek · April 30, 2020, 9:21am

Hi,

There is really no such thing as “CVX ignores constraint” or “Solver ignores constraint”. They never do (unless it is some horrible bug). In almost all cases a constraint is violated for numerical reasons or there is a bug in your model.

Therefore it is necessary in such cases to see solver log output, or have a fully reproducible example with input data or anything that could help evaluate what sort of numerical trouble your problem gets into. Is the constraint violation substantial or is it just numerical noise? Have you tried various solvers?

Mobin · April 30, 2020, 9:40am

Hi Michal.

Thanks for fast response. The violation is really large, e.g. 5000 %. SDP3 can solve the problem but in an infeasible path-following algorithms. Sedumi cannot solve it and the status is infeasible. These are for both versions 2.1 and 2.2.

The code is ready to be run except two data packs, which I can send to you if you give me your mail address, please.

Michal_Adamaszek · April 30, 2020, 10:06am

Please post the log outputs. That should clarify some things.

If one solver returns a bad solution and the other says infeasible, isn’t that already a huge enough warning sign? Try also Mosek if you can.

Mobin · April 30, 2020, 10:25am

Sorry, I don’t have the licence for commercial solvers.

It is the log output without any warning:

Calling SDPT3 4.0: 21383 variables, 719 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 719
 dim. of sdp    var  = 888,   num. of sdp  blk  = 96
 dim. of linear var  = 8969
 dim. of free   var  = 294
 566 linear variables from unrestricted variable.
 *** convert ublk to lblk
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|9.0e+04|9.3e+01|1.2e+12| 3.998214e+08  0.000000e+00| 0:0:00| spchol 
  SMW too ill-conditioned, switch to LU factor, Inf.
  warning: symqmr failed: 0.3 
  switch to LU factor. splu 294  294 
  linsysolve: solution contains NaN or inf
  stop: difficulty in computing corrector directions
-------------------------------------------------------------------
 number of iterations   =  1
 primal objective value =  3.99821360e+08
 dual   objective value =  0.00000000e+00
 gap := trace(XZ)       = 1.17e+12
 relative gap           = 2.94e+03
 actual relative gap    = 1.00e+00
 rel. primal infeas (scaled problem)   = 9.04e+04
 rel. dual     "        "       "      = 9.29e+01
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 8.9e+03, 0.0e+00, 1.3e+08
 norm(A), norm(b), norm(C) = 3.0e+03, 1.4e+00, 1.4e+06
 Total CPU time (secs)  = 1.19  
 CPU time per iteration = 1.19  
 termination code       =  0
 DIMACS: 1.2e+05  0.0e+00  2.5e+03  0.0e+00  1.0e+00  2.9e+03
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +128.174

Erling · April 30, 2020, 10:29am

There is a warning.

0|0.000|0.000|9.0e+04|9.3e+01|1.2e+12| 3.998214e+08 0.000000e+00| 0:0:00| spchol
SMW too ill-conditioned, switch to LU factor, Inf.
warning: symqmr failed: 0.3

So you might have linear dependencies in your constraints since it happens in the first iteration.

Michal_Adamaszek · April 30, 2020, 10:36am

As you can see the solver gives up immediately and returns completely bad solution. In almost every line it tells you something went wrong. Probably it is a garbage-in garbage-out situation. You need to work on the scaling and other numerics, check your data, maybe solve smaller instances to try things out.

It is “well known” in some circles that certain power flow datasets out there require a lot of cleaning before they can be fed in a solver. Maybe your case is similar.

Mobin · April 30, 2020, 11:18am

Thank you both. Something came to my mind about the input data and its constraining. I will try to solve the problem. Probably, I bother you again .

Mark_L_Stone · April 30, 2020, 12:16pm

It is interesting (concerning) that the solver seems to end in a horrible way, yet CVX reports Solved and produces an Optimal value. It seems to me that CVX should have produced Status Failed, and populated all variables with NaN.

CVX 3.0beta is full of bugs and should not be used.