%% EXACT CONVEX RELAXATION OF OPTIMAL POWER FLOW IN RADIAL NETWORKS
% Submitted by : Prajakta Muley and Mehroze Jalal
%% PROBLEM FORMULATION
% The optimal power flow (OPF) problem aims to minimize the power loss and generation cost and hence the system cost of a radial distribution network. The approach used to deal with 
% this challenge convexifies the constraints imposed by Kirchoff's laws. 
%  
% For radial networks three types of sufficient conditions have been developed that guarantee the exactness of convex relaxations: 
%
% 1. Power Injections, 
% 2. Voltage Angles,
% 3. Voltage Magnitudes
%
% The main goal is to prove that after modifying the OPF problem for radial networks slightly, the corresponding SOCP relaxation is exact under 
% a mild condition that can be checked a priori. 

%% PROPOSED SOLUTION 
% To solve the above stated problem, this code demonstrates the solutions to the following optimization problems: 
%  
% 1. The original OPF problem
%
% <<PROJECT11.PNG>>
% 
% <<PROJECT1_2.PNG>>
%  
% The voltage deviation is set to 5% from their nominal values, hence 0.95^2 <= vi <=1.05^2. The function fi denotes the generation cost at bus i. The paper states that if fi (x) = x, 
% then the objective function is the total power loss on the network. The OPF problem minimizes the generation cost subject to power flow constraints, power injection constraints 
% and voltage constraints.  
% 
% 2. The original SOCP problem
%
% The third equality constraint above is nonconvex and can be relaxed to inequality constraint to obtain the following second order cone program.  SOCP is exact if every of its optimal solutions 
% satisfies the third constraint in the OPF problem. 
% In this model, we cosider the same constraints as the OPF problem but only one constraint will change which is mentioned below:
%
% <<PROJECT1_3.PNG>>
% 
% s.t. 
% 
% same 1, 2, 4 and 5 constraints and 
%
% <<PROJECT1_4.PNG>>
%
% 3. The modified OPF problem 
%
% Using two constraints, C1 and C2 the OPF problem on network distribution systems can be modified slightly. The constraints are shown below:  
%
%
% <<PROJECT1_5.PNG>>
% 
% C1 can be defined from the network model which is given below:
% 
% <<PROJECT1_7.PNG>>
%
% svolt can be defined as
%
% 
% <<PROJECT1_6.PNG>>
%
% Using the Linear DistFlow Model, and its solutions the Svolt constraint (C2) above is defined. 
% 
% <<PROJECT1_8.PNG>>
%
% Hence the modified OPF problem is : 
% 
% <<PROJECT1_9.PNG>>
% 
% Here also we consider same constraints as normal OPF problem.
% 
% <<PROJECT1_10.PNG>>
% 
% The modified SOCP problem 
%
% A modification to OPF is necessary to ensure an exact SOCP, since otherwise examples exist where SOCP is not exact. The feasible sets of OPF-m and OPF are similar since both the vi ranges used 
% in each are close to each other in practice. 
% 
% <<PROJECT1_11.PNG>>
%
% Since OPF-m is obtained by imposing additional constraints on OPF, C2. It follows that SOCP-m relaxation is exact under C1, a mild condition that can be checked a priori. 
%
%% DATA SOURCE
% The IEEE 14-bus network data has been utilized.
% 
POWERFLOWS= [
%Bus         V           phase       P gen       Q gen       P load      Q load                 
  01      1.06        0           2.3258     -0.14978     0           0          
  02      1.045      -0.0871      0.4         0.48824     0.217       0.127      
  03      1.01       -0.22267     0           0.27372     0.942       0.19       
  04      1.012      -0.1785      0           0           0.478       0.04                 
  05      1.016      -0.15273     0           0           0.076       0.016      
  06      1.07       -0.25161     0           0.2251      0.112       0.075      
  07      1.06       -0.245       0.54406    -0.1454      0.54406    -0.1454
  08      1.09       -0.23091     0           0.25163     0           0          
  09      1.0328     -0.25853     0           0           0.295       0.166      
  10      1.0318     -0.26223     0           0           0.09        0.058      
  11      1.0471     -0.25897     0           0           0.035       0.018      
  12      1.0534     -0.26645     0           0           0.061       0.016      
  13      1.047      -0.2671      0           0           0.135       0.058      
  14      1.0207     -0.28018     0           0           0.149       0.05];

%         |  From |  To   |   R     |   X     |     B/2  |  X'     |
%         |  Bus  | Bus   |  pu     |  pu     |     pu   | TAP (a) |
LINEDATA =  [1      2       0.01938   0.05917    0.0264         1
             1      5       0.05403   0.22304    0.0246         1
             2      3       0.04699   0.19797    0.0219         1
             2      4       0.05811   0.17632    0.0170         1
             3      4       0.06701   0.17103    0.0064         1
             4      5       0.01335   0.04211    0.0            1
             6     11       0.09498   0.19890    0.0            1
             6     12       0.12291   0.25581    0.0            1
             6     13       0.06615   0.13027    0.0            1
             9     10       0.03181   0.08450    0.0            1
             9     14       0.12711   0.27038    0.0            1
            10     11       0.08205   0.19207    0.0            1
            12     13       0.22092   0.19988    0.0            1
            13     14       0.17093   0.34802    0.0            1 ];
        
% we solved this model in PSAT for power flow problem
% From PSAT(power system analysis toolbox), we get the P flow and Q flow values, 

% frombus  tobus      Line        P Flow      Q Flow      P Loss      Q Loss
 LINEFLOWS=    [01      02      10          1.5712     -0.2046      0.04311     0.07312
                01      05      13          0.7546      0.05482     0.02771     0.06135    
                02      03      11          0.73462     0.03538     0.02337     0.05222                                    
                02      04      15          0.55939     0.01594     0.01672     0.01117  
                03      04      12         -0.23076     0.06689     0.00397    -0.02524
                04      05      14         -0.59712     0.09398     0.00478     0.00192
                06      11      5           0.0818      0.08439     0.00115     0.0024
                06      12      2           0.08037     0.03119     0.0008      0.00166        
                06      13      4           0.18272     0.09743     0.00248     0.00488       
                09      10      7           0.04487    -0.00475     6e-05       0.00016
                09      14      8           0.08719     0.00597     0.00091     0.00194
                10      11      6          -0.04519    -0.06291     0.00046     0.00108
                12      13      3           0.01857     0.01353     0.00011     0.0001 
                13      14      9           0.06371     0.04799     0.00099     0.00202];

% Computational time SDP =  0.138s
%% NUMERICAL SOLUTION 
%
% Two tests were performed to demonstrate that: 
%
% 1. SOCP is simpler computationally than SDP and  
%
% 2. The feasible sets of OPF and OPF-m are similar. We define modification gap that quantifies the difference between the feasible sets of OPF  and OPF-m, and show that this gap is small. 
%
% For the first test, the computation time for SDP was obtained from the reference  and the one for SOCP was computed using the MATLAB code.
%
% The second test shows that OPF-m eliminates some feasible points of OPF that are close to the voltage upper bounds for all test networks. 
% 
% 
% <<PROJECT1_12.PNG>>
% 
% <<PROJECT1_13.PNG>>
%  
% <<PROJECT1_14.PNG>>
% 
% <<PROJECT1_15.PNG>>
%
% If error is small, then FOPF-m is similar toFOPF. Any point w that is feasible for OPF but infeasible for OPF-m is close to the voltage upper bound since vi > vi ? error for some i which is belongs to Nplus set. Such points are perhaps undesirable for robust operation.
%% CODE AND DEMONSTRATION OF RESULTS

Pij=LINEFLOWS(:,4); % active power flow from bus i to j
Qij=LINEFLOWS(:,5); % reactive power flow from bus i and j
Sij=Pij+1i*Qij; % sending end power flow from bus i and j 
Sijmag= sqrt((Pij).^2+(Qij).^2); % magnitude of Sij


vi=POWERFLOWS(:,2); % magnitude of complex voltage sqaured

for n=1:13
    vj=vi(n+1);
end


rij=LINEDATA(:,3); % resitance of line between bus i and j
xij = LINEDATA(:,4);% reactance of line between bus i and j
zij=rij+1i*xij; % impedance of line between bus i and j
zijmag=sqrt((rij).^2+(xij).^2); % magnitude of zij

pi=POWERFLOWS(:,6); % real power injection of bus i
qi=POWERFLOWS(:,7); % reactive power ijection of bus i
si=pi+1i*qi; % total power injection of bus i
gij=1./rij; % condusctane of line between bus i and j

lij=((Pij).^2+(Qij).^2)/vi; % square of magnitude of complex current Iij
lij2 = lij(:,8); % square of magnitude of complex current Iij

h=[11]; % downstream bus (given that h belongs to set of buses i)
%hi defines from bus h to bus i
Phi= [0 0 0 0 0 0.0818 0 0 0 -0.04519 0 0 0 0]'; 
Qhi= [0 0 0 0 0 0.08439 0 0 0 -0.06291 0 0 0 0]';
Shi=Phi+1i*Qhi;
rhi=[0 0 0 0 0 0.09498 0 0 0 0.08205 0 0 0 0]';
xhi=[0 0 0 0 0 0.19890 0 0 0 0.19207 0 0 0 0]';
zhi=rhi+1i*xhi;

lhi= abs(Shi).^2/vi;
lhi2=lhi(:,8);


x=Shi-(zhi.*lhi2);

ph=[0 0 0 0 0 0 0 0 0 0 0.035 0 0 0]';
qh=[0 0 0 0 0 0 0 0 0 0 0.018 0 0 0]';
sh = ph+ 1i*qh;
Sijcap= sh;
% bus k is upstream bus such that it belongs to set of bus i
% jk = from bus j to bus k
Sjkcap=Sijcap;
zjk=conj(zij);
rjk=real(zij);
Pjkcap=real(Sjkcap);
% vo is substation bus which means reference bus/ slack bus
vo=1;
vicap=vo+2.*(rjk.*Pjkcap);
%ORIGINAL OPF
tic
cvx_begin
cvx_solver sedumi

variable f1(14)

minimize sum(f1)
subject to

    f1==pi;
    Sij==si+x;
    vi- vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2==Sijmag./vi;
    0.9025<=vi<=1.1025;
    
    
cvx_end
cvx_optval
toc

%ORIGINAL SOCP
tic
cvx_begin
cvx_solver sedumi

variable f2(14)

minimize sum(f2)
subject to

    f2==pi;
    Sij==si+x;
    vi- vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2>=Sijmag./vi;
    si == -pi - 1i*qi;
    0.81<= vi<= 1.21;
    
cvx_end
cvx_optval
toc

%OPF-m

tic
cvx_begin
cvx_solver sedumi

variable f3(14)

minimize sum(f)
subject to

    f3== pi;
    Sij==si+x;
    vi- vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2==Sijmag./vi;
    0.8<=vicap<=1.21;
    0.8<=vi<=1.21;
    
cvx_end
cvx_optval
toc

%SOCP-m
tic
cvx_begin
cvx_solver sedumi

variable f4(14)

minimize sum(f4)
subject to

    f4== pi;
    Sij==si+x;
    vi-vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2>=Sijmag./vi;
    0.81<=vi<=1.21;
    
cvx_end
cvx_optval
toc

%TEST2
cvx_begin
cvx_solver sedumi


E=norm(vicap-vi,Inf);


maximize E
subject to 
    Sij==si+x;
    vi- vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2==Sijmag./vi;
    0.9025<=vi<=1.1025;
    
cvx_end
cvx_optval 


cvx_begin
cvx_solver sedumi

variable f5(14)

minimize sum(f5)
subject to

    f5== pi;
    Sij==si+x;
    vi- vj==2.*(rij.*Pij)- ((zijmag).^2).*lij2;
    lij2==Sijmag./vi;
    0.81<=vi<=1.21-E;
cvx_end
cvx_optval

%% SUMMARY AND CONCLUSION
% We have proved that SOCP is exact if condition C1 and C2 hold. We have proposed a modified OPF problem that includes the additional constraint that power injections lie in Svolt, such that the corresponding SOCP relaxation is exact under C1. 
% Empirical studies show that the feasible sets (vicap<=vi) of OPF and OPF-m are close for the IEEE 14 bus network. We also proved that SOCP has at least one optimal solution if it is convex and convex. 
 


