%% Optimal Inverter VAR Control in Distribution Systems with High PV Penetration
%
%% Problem Formulation
%
% The aim of the study is to demonstrate the benefits of inverter reactive
% power control to mitigate the rapid and large voltage fluctuations due to
% the high penetration of photovoltaic generation and the resulting reverse
% power flow.
%

%% Proposed Solution
%
%
% The approach for this problem is formulated as a radial optimal power
% flow (OPF) to minimiza the loss and energy consumptions, subject to
% constraints on voltage magnitudes.
%
% A minimization of power loss is formulated as in the paper page 4, equation (17), which considers 
% the line losses as the first part, the Conservation Voltage Reduction
% effects (CVR) as the second part, and the losses in the inverter due to
% reactive power fow as the thord part.
%
% The minimization problem is subjected to a bunch of voltage and power
% contraint in the while system, which consists of all the power, current,
% and voltage equations, which are shown in the paper page 4, equation (18)-(27).
%
% Tn ths study, first, we investigate the interaction between the solar
% plant reactive power, active power, and overall load under the voltage
% constraint.
%
% In this algorithm operation, we first intialize all the system
% parameters, then real-life data from South California Edison (SCE) company are
% applied to do the simulation. The optimization algorithm is solved by
% the CVX solver.

%% System parameter initialization
%
% First, all the parameters will be initiallized to provide a general
% picture of the problem.
%

pcs=[0.2 0.3 0.4];              % loop value for system load
pgs_45 = [1:0.5:4];           % input parameter of PV real power injection
qgs_45 = zeros(size(pgs_45));

for pc_s=1:3
    
for pg_i = 1:7



i_total=56;                 % define system size
j_total=56;     

pc=zeros(i_total,1);         % real power demand
qc=zeros(i_total,1);         % reactive power demand
sc=zeros(i_total,1);         % total power demand
qsc=zeros(i_total,1);        % shunt capacitor vector
pg=zeros(i_total,1);         % real power generated by PVs and inverters
qg=zeros(i_total,1);         % reactive power generated by PVs and inverters

r=zeros(i_total,j_total);    % the line resistance
x=zeros(i_total,j_total);    % the line reactance

l=zeros(i_total,j_total);    % current relaxation
v=12*ones(i_total,1);        % voltage relaxation
s=zeros(i_total,1);          % power relaxation
t=zeros(i_total,1);          

alpha=zeros(i_total,1);      % CVR coefficient
cv=ones(i_total,1);          % inverter loss coefficients
cr=ones(i_total,1);         
qg_max=zeros(i_total,1);     % optimization constraints parameters
V_min=zeros(i_total,1);
V_max=zeros(i_total,1);
pc_min=zeros(i_total,1);
qc_min=zeros(i_total,1);

P=zeros(i_total,j_total);    %real power flow 
Q=zeros(i_total,j_total);    %reactive power flow
E=zeros(i_total,j_total);    %network connection (graph arc)


%% Data import of the simulation
% In this study, we input the historical SCADA data from SCE company into a
% chosen SCE's 56-node distribution feeders with very high penetration of
% photovoltaics.

E(1,2)=1;                   % E Matrix, imply the connection between each node, 1:connected, 0: open
E(2,3)=1;
E(2,4)=1;
E(4,5:7)=1;
E(7,8)=1;
E(8,9:10)=1;
E(10,11)=1;
E(11,12:13)=1;
E(13,14:15)=1;
E(15,16:17)=1;
E(17,18)=1;
E(18,19)=1;
E(4,20)=1;
E(20,21)=1;
E(21,22)=1;
E(20,23)=1;
E(23,24:25)=1;
E(25,26)=1;
E(26,27)=1;
E(27,28)=1;
E(28,29)=1;
E(29,30)=1;
E(30,31)=1;
E(26,32)=1;
E(32,33)=1;
E(32,34)=1;
E(34,35:36)=1;
E(36,37)=1;
E(34,38)=1;
E(38,39)=1;
E(39,40)=1;
E(34,41)=1;
E(41,42)=1;
E(42,43:46)=1;
E(41,47)=1;
E(47,48:49)=1;
E(49,50:51)=1;
E(51,52:53)=1;
E(53,54:56)=1;


r(1,2)=0.16;             % r matrix, imply network resistances
r(2,3)=0.824;
r(2,4)=0.144;
r(4,5)=1.026;
r(4,6)=0.741;
r(4,7)=0.528;
r(7,8)=0.358;
r(8,9)=2.032;
r(8,10)=0.502;
r(10,11)=0.372;
r(11,12)=1.431;
r(11,13)=0.429;
r(13,14)=0.671;
r(13,15)=0.457;
r(15,16)=1.008;
r(15,17)=0.153;
r(17,18)=0.971;
r(18,19)=1.885;
r(4,20)=0.138;
r(20,21)=0.251;
r(21,22)=1.818;
r(20,23)=0.225;
r(23,24)=0.127;
r(23,25)=0.284;
r(25,26)=0.171;
r(26,27)=0.414;
r(27,28)=0.21;
r(28,29)=0.395;
r(29,30)=0.248;
r(30,31)=0.279;
r(26,32)=0.205;
r(32,33)=0.263;
r(32,34)=0.071;
r(34,35)=0.625;
r(34,36)=0.51;
r(36,37)=2.018;
r(34,38)=1.062;
r(38,39)=0.61;
r(39,40)=2.349;
r(34,41)=0.115;
r(41,42)=0.159;
r(42,43)=0.934;
r(42,44)=0.506;
r(42,45)=0.095;
r(42,46)=1.915;
r(41,47)=0.157;
r(47,48)=1.641;
r(47,49)=0.081;
r(49,50)=1.727;
r(49,51)=0.112;
r(51,52)=0.674;
r(51,53)=0.07;
r(53,54)=2.041;
r(53,55)=0.813;
r(53,56)=0.141;


x(1,2)=0.388;           % x matrix, imply network reactance
x(2,3)=0.315;
x(2,4)=0.349;
x(4,5)=0.421;
x(4,6)=0.466;
x(4,7)=0.468;
x(7,8)=0.314;
x(8,9)=0.798;
x(8,10)=0.441;
x(10,11)=0.327;
x(11,12)=0.999;
x(11,13)=0.377;
x(13,14)=0.257;
x(13,15)=0.401;
x(15,16)=0.385;
x(15,17)=0.134;
x(17,18)=0.722;
x(18,19)=0.721;
x(4,20)=0.334;
x(20,21)=0.096;
x(21,22)=0.695;
x(20,23)=0.542;
x(23,24)=0.028;
x(23,25)=0.687;
x(25,26)=0.414;
x(26,27)=0.386;
x(27,28)=0.196;
x(28,29)=0.369;
x(29,30)=0.232;
x(30,31)=0.26;
x(26,32)=0.495;
x(32,33)=0.073;
x(32,34)=0.171;
x(34,35)=0.273;
x(34,36)=0.209;
x(36,37)=0.829;
x(34,38)=0.406;
x(38,39)=0.238;
x(39,40)=0.964;
x(34,41)=0.278;
x(41,42)=0.384;
x(42,43)=0.383;
x(42,44)=0.163;
x(42,45)=0.195;
x(42,46)=0.769;
x(41,47)=0.379;
x(47,48)=0.67;
x(47,49)=0.196;
x(49,50)=0.709;
x(49,51)=0.27;
x(51,52)=0.275;
x(51,53)=0.17;
x(53,54)=0.78;
x(53,55)=0.334;
x(53,56)=0.34;

% load matrix
pf = 0.9;                   % power factor
pfsin = sqrt(1-pf^2);

sc(1)=-4;
sc(3)=0.057;
sc(5)=0.121;
sc(6)=0.049;
sc(7)=0.053;
sc(8)=0.047;
sc(9)=0.068;
sc(10)=0.048;
sc(11)=0.067;
sc(12)=0.094;
sc(14)=0.057;
sc(16)=0.053;
sc(17)=0.057;
sc(18)=0.112;
sc(19)=0.087;
sc(22)=0.063;
sc(24)=0.135;
sc(25)=0.1;
sc(27)=0.48;
sc(28)=0.38;
sc(29)=0.044;
sc(31)=0.053;
sc(32)=0.223;
sc(33)=0.123;
sc(34)=0.067;
sc(35)=0.094;
sc(36)=0.097;
sc(37)=0.281;
sc(38)=0.117;
sc(39)=0.131;
sc(40)=0.03;
sc(41)=0.046;
sc(42)=0.054;
sc(43)=0.083;
sc(44)=0.057;
sc(46)=0.134;
sc(47)=0.045;
sc(48)=0.196;
sc(50)=0.045;
sc(52)=0.315;
sc(54)=0.061;
sc(55)=0.055;
sc(56)=0.13;

qg_max(45)=5;                   % reactive power upper bound

for i=1:i_total                 % calculate the real and reactive power
    pc(i) = pcs(pc_s)*sc(i)*pf;
    qc(i) = pcs(pc_s)*sc(i)*pfsin;
end

for i=1:i_total
    alpha(i)=pc(i);
end

qg(45)=0;                       % PV reactive power initialization


    pg(45) = pgs_45(pg_i);      % PV real power initialization
    
%% Solution
%
% The optimal power flow problem is sovled by CVX optimization toolbox.


i=0;
k=0;

cvx_begin                       % start the optimization tools

variable l(i_total,j_total);    % state all the optimization variables
variable v(i_total,1);
variable s(i_total,1);       
variable t(i_total,1);
variable P(i_total,j_total);
variable Q(i_total,j_total);
variable qg(i_total,1);


minimize (trace(r'*l)+alpha'*v+(cv'*s+cr'*t))       % the total minimization objective, eqa(17) in paper page 4

subject to


    s >= norm([pg qg]);                             % s constraint, eqa(22)
    t >= pg.^2+qg.^2;                               % t constraint, eqa(23) 
    abs(qg) <= qg_max;                              % PV reactive power constraint, eqa(24)
    4 <= v <= 18;                                   % V upper and down constraint, eqa(25)
    for h=2:i_total
        pc_min <= pc(h);                            % demand power constraint, eqa(26)
        qc_min <= qc(h);
    end
    
    for i=1:i_total
        for k=i+1:j_total
            P(i,k) == P(k,:)*E(k,:)' + r(i,k)*l(i,k)+ pc(k)-pg(k);                         % P constraint, eqa(18)
            Q(i,k) == Q(k,:)*E(k,:)' + x(i,k)*l(i,k)+ qc(k)-qg(k)-v(k)*qc(k);              % Q constraint, eqa(19)
            
            v(k)*E(i,k) == (v(i)-2*(r(i,k)*P(i,k)+x(i,k)*Q(i,k))+(r(i,k)^2+x(i,k)^2)*l(i,k))*E(i,k);        % V constraint, eqa(20)
            norm([2*P(i,k); 2*Q(i,k); l(i,k)-v(i)])*E(i,k) <= (l(i,k)+v(i))*E(i,k);                         % lij constraint, eqa(21)
        end
    end

    
cvx_end                         % end the optimization tools
    
    qgs_45(pg_i) = qg(45);
end
%% Visulization of results
%
% In this section, we show the plot of the relation between photovaltics
% output real power and reactive power (VAR injection).
% The first plot shows the relation when systme load is 20% of the peak
% load, the second one shows 30% of the peak load, and the third one shows
% 40%.

    figure;
    plot(pgs_45,qgs_45,'bo-');
    grid on;
    xlabel('solar real power [MW]');
    ylabel('Optimal inverter VAR injection [MW]');
    title('Optimal inverter reactive power vs power output');

end


%% Analysis and Conclusions
%
% By choosing appropriate parameters,we obtained results similar to the
% original paper. 
% However, there are some data which are not announced clearly in paper. In
% order to make it work, we made some changes, shown as follow:
%
% 1) The parameter cv and cr are not given in the paper, we simply set
% cv=cr=1;
%
% 2) The first node(#1)is supposed to be feed from a susbstation,by means
% of which it's able to maintain a constant voltage value. Thus the load 
% demand of node #1 should be minus. However there is nothing mentioned
% about the supply in this article.
%
% 3) Due to the limitted data access of the simulation, we loose the
% voltage constraints a little to make CVX feasible.
%
% 4)The plot from the visulization part fits the Fig.4 in paper. As for
% Fig.5, since we cannot access to the value setting of cv and cr, the
% simulation will not be necessary.
%
% 5) Due to the limitation of time, we didn't make the plot more smooth, as
% the current point number already cost too much times for one
% simulation(almost two hours).