%% Minimum Spanning Tree Calculation for the 14-Node system of [1].
% For a graphic view of the problem, look at the accompanying image file.
% This work has been submitted for publication to IET Generation,
% Transmission & Distribution, Jan. 2014.
% You need CVX package to run the file. It is freely available online at:
% http://cvxr.com/cvx/download/
% [1] S. Civanlar, J. J. Grainger, H. Yin, and S. S. H. Lee, “Distribution 
%     feeder reconfiguration for loss reduction,” IEEE Trans. Power Del.,
%     vol.3, no.3, pp. 1217–1223, Jul. 1988.
clc
clear
% Per-unit values of branch current magnitues:
Iij = [
%  From To    |I| 
    1	2	0.110828547
    2	3	0.059656128
    2	4	0.028396877
    4	5	0.012727385
    1	6	0.112127193
    6	7	0.074301973
    6	8	0.006425698
    7	9	0.023429612
    7	10	0.048638608
    1	11	0.076275149
    11	12	0.025416302
    11	13	0.039530604
    13	14	0.029821935
    3	9	0.028983145
    8	12	0.020056312
    5	14	0.008132599    
];
% Numbe of nodes
n = max(max(Iij(:,1:2)));
% Number of branches
m = length(Iij(:,1));
% Edge Weights
w = sparse(Iij(:,1),Iij(:,2),Iij(:,3));w(n,:) = 0;
w = full(w);
w = w'+w;
% Beging formulating Linear Programming problem in CVX
cvx_begin
% cvx_solver SDPT3
cvx_solver sedumi
cvx_precision high
variable x(n,n)
variables yab yba yac yca ybc ycb yad_1 yad_2 yad_3 yad_4 yad_5 ybd_1 ybd_2 ybd_3 ybd_4 ycd_1 ycd_2 ycd_3
maximize(sum(sum(w.*x)))
subject to
for i=1:n-1
    sum(x(i,w(:,i)>0)) == 1; %#ok<*EQEFF>
end
yab+yac+yad_1+yad_2+yad_3+yad_4+yad_5 == 1;
yba+ybc+ybd_1+ybd_2+ybd_3+ybd_4 == 1;
yca+ycb+ycd_1+ycd_2+ycd_3 == 1;
x(1,2)+x(2,1)+yab+yba == 1;
x(2,3)+x(3,2)+ybd_1 == 1;
x(2,4)+x(4,2)+yad_1 == 1;
x(4,5)+x(5,4)+yad_2 == 1;
x(1,6)+x(6,1)+ybc+ycb == 1;
x(6,7)+x(7,6)+ybd_4 == 1;
x(6,8)+x(8,6)+ycd_3 == 1;
x(7,9)+x(9,7)+ybd_3 == 1;
x(7,10)+x(10,7) == 1;
x(1,11)+x(11,1)+yac+yca == 1;
x(11,12)+x(12,11)+ycd_1 == 1;
x(11,13)+x(13,11)+yad_5 == 1;
x(13,14)+x(14,13)+yad_4 == 1;
x(3,9)+x(9,3)+ybd_2 == 1;
x(8,12)+x(12,8)+ycd_2 == 1;
x(5,14)+x(14,5)+yad_3 == 1;
x(w>0) <= 1;x(w>0) >= 0;x(w==0) == 0; %#ok<*VUNUS>
yab>=0;yba>=0;yac>=0;yca>=0;ybc>=0;ycb>=0;yad_1>=0;yad_2>=0;yad_3>=0;yad_4>=0;yad_5>=0;ybd_1>=0;ybd_2>=0;ybd_3>=0;ybd_4>=0;ycd_1>=0;ycd_2>=0;ycd_3>=0;
yab<=1;yba<=1;yac<=1;yca<=1;ybc<=1;ycb<=1;yad_1<=1;yad_2<=1;yad_3<=1;yad_4<=1;yad_5<=1;ybd_1<=1;ybd_2<=1;ybd_3<=1;ybd_4<=1;ycd_1<=1;ycd_2<=1;ycd_3<=1;
cvx_end
% End of the optimization routine
% Calculating the values for x
x = double(round(x));
% Defining connected branches in the minimum spanning tree
u = zeros(n,n);
for i = 1:n
    for j = 1:n
        u(i,j) = x(i,j)+x(j,i);
    end
end
% Defining disconnected branches in the minimum spanning tree
u_open = zeros(n,n);
for i = 1:n
    for j = 1:n
        if u(i,j)<0.5 && w(i,j)>0
            u_open(i,j) = 1;
        end
    end
end
u_open = sparse(u_open);


