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:TM1
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(i4);
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_d1P(8,:);
P_ord2 = P_d2P(9,:);
P_ord3 = P_d3P(10,:);
P_critic = P_crP(11,:);
P_pv1 = Ps/3P(16,:);
P_pv2 = Ps/3P(57,:);
P_pv3 = Ps/3P(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/3P(17,:); % Produced active power of the wind generator
P_wind2 = Pw/3P(59,:);
P_wind3 = Pw/3P(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(i18,t)<=P_gmax;
0<=P(i,t)<=Q_gmax;
(P(i,t)^2)+(P(i18,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_ord3P_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_d3Q_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_11Y1_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)]; % Threephase 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