%% Minimum Spanning Tree Calculation for the 33-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] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution 
%     systems for loss reduction and load balancing,” IEEE Trans. Power 
%     Del., vol.4, no.2, pp. 1401–1407, Apr. 1989.
clc
clear
% Per-unit values of branch current magnitues:
Iij = [
%  From To    |I| 
    1	2	4.520470523
    2	3	3.150129873
    3	4	1.439156927
    4	5	1.293102337
    5	6	1.224601742
    6	7	0.469050527
    7	8	0.249060881
    8	9	0.420748772
    9	10	0.073555047
    10	11	0.032591987
    11	12	0.07153755
    12	13	0.291691526
    13	14	0.220454879
    14	15	0.074147536
    15	16	0.319270551
    16	17	0.263243697
    17	18	0.212827232
    2	19	1.256974466
    19	20	1.160217422
    20	21	1.062400281
    21	22	0.527598591
    3	23	1.614266709
    23	24	1.50954017
    24	25	1.036189266
    6	26	0.714211685
    26	27	0.649602599
    27	28	0.585455102
    28	29	0.525255147
    29	30	0.955348566
    30	31	0.448455374
    31	32	0.285629473
    32	33	0.126188837
    8	21	0.437493453
    9	15	0.300542463
    12	22	0.429803323
    18	33	0.149966278
    25	29	0.569815961
];
% 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 y1_6_1 y1_6_2 y1_6_3 y1_6_4 y1_4_1 y1_4_2 y1_4_3 y1_5_1 y1_5_2 y1_2
variables y2_1 y2_6_1 y2_6_2 y2_3_1 y2_3_2 y2_3_3 y2_5
variables y3_5 y3_2_1 y3_2_2 y3_2_3 y3_6_1 y3_6_2 y3_6_3
variables y4_6_1 y4_6_2 y4_6_3 y4_6_4 y4_5_1 y4_5_2 y4_5_3 y4_5_4 y4_1_1 y4_1_2 y4_1_3
variables y5_6_1 y5_6_2 y5_6_3 y5_6_4 y5_6_5 y5_6_6 y5_6_7 y5_6_8 y5_4_1 y5_4_2 y5_4_3 y5_4_4 y5_1_1 y5_1_2 y5_2 y5_3
maximize(sum(sum(w.*x)))
subject to
for i=1:n-1
    sum(x(i,w(:,i)>0)) == 1; %#ok<*EQEFF>
end

y1_6_1+y1_6_2+y1_6_3+y1_6_4+y1_4_1+y1_4_2+y1_4_3+y1_5_1+y1_5_2+y1_2 == 1;
y2_6_1+y2_6_2+y2_3_1+y2_3_2+y2_3_3+y2_1+y2_5 == 1;
y3_6_1+y3_6_2+y3_6_3+y3_2_1+y3_2_2+y3_2_3+y3_5 == 1;
y4_6_1+y4_6_2+y4_6_3+y4_6_4+y4_1_1+y4_1_2+y4_1_3+y4_5_1+y4_5_2+y4_5_3+y4_5_4 == 1;
y5_2+y5_3+y5_6_1+y5_6_2+y5_6_3+y5_6_4+y5_6_5+y5_6_6+y5_6_7+y5_6_8+y5_4_1+y5_4_2+y5_4_3+y5_4_4+y5_1_1+y5_1_2 == 1;

x(1,2)  +x(2,1)                    == 1;
x(2,3)  +x(3,2)    +y1_6_1         == 1;
x(3,4)  +x(4,3)    +y1_4_1+y4_1_1  == 1;
x(4,5)  +x(5,4)    +y1_4_2+y4_1_2  == 1;
x(5,6)  +x(6,5)    +y1_4_3+y4_1_3  == 1;
x(6,7)  +x(7,6)    +y1_5_1+y5_1_1  == 1;
x(7,8)  +x(8,7)    +y1_5_2+y5_1_2  == 1;
x(8,9)  +x(9,8)    +y2_5+y5_2      == 1;
x(9,10) +x(10,9)   +y2_3_1+y3_2_1  == 1;
x(10,11)+x(11,10)  +y2_3_2+y3_2_2  == 1;
x(11,12)+x(12,11)  +y2_3_3+y3_2_3  == 1;
x(12,13)+x(13,12)  +y3_6_1         == 1;
x(13,14)+x(14,13)  +y3_6_2         == 1;
x(14,15)+x(15,14)  +y3_6_3         == 1;
x(15,16)+x(16,15)  +y5_6_8         == 1;
x(16,17)+x(17,16)  +y5_6_7         == 1;
x(17,18)+x(18,17)  +y5_6_6         == 1;
x(2,19) +x(19,2)   +y1_6_2         == 1;
x(19,20)+x(20,19)  +y1_6_3         == 1;
x(20,21)+x(21,20)  +y1_6_4         == 1;
x(21,22)+x(22,21)  +y2_6_1         == 1;
x(3,23) +x(23,3)   +y4_6_1         == 1;
x(23,24)+x(24,23)  +y4_6_2         == 1;
x(24,25)+x(25,24)  +y4_6_3         == 1;
x(6,26) +x(26,6)   +y5_4_1+y4_5_1  == 1;
x(26,27)+x(27,26)  +y5_4_2+y4_5_2  == 1;
x(27,28)+x(28,27)  +y5_4_3+y4_5_3  == 1;
x(28,29)+x(29,28)  +y5_4_4+y4_5_4  == 1;
x(29,30)+x(30,29)  +y5_6_1         == 1;
x(30,31)+x(31,30)  +y5_6_2         == 1;
x(31,32)+x(32,31)  +y5_6_3         == 1;
x(32,33)+x(33,32)  +y5_6_4         == 1;
x(21,8) +x(8,21)   +y1_2+y2_1      == 1;
x(9,15) +x(15,9)   +y3_5+y5_3      == 1;
x(22,12)+x(12,22)  +y2_6_2         == 1;
x(33,18)+x(18,33)  +y5_6_5         == 1;
x(25,29)+x(29,25)  +y4_6_4         == 1;

x(w>0) <= 1;x(w>0) >= 0;x(w==0) == 0; %#ok<*VUNUS>

y1_6_1>=0;y1_6_2>=0;y1_6_3>=0;y1_6_4>=0;y1_4_1>=0;y1_4_2>=0;y1_4_3>=0;y1_5_1>=0;y1_5_2>=0;y1_2>=0;
y2_1>=0;y2_6_1>=0;y2_6_2>=0;y2_3_1>=0;y2_3_2>=0;y2_3_3>=0;y2_5>=0;
y3_5>=0;y3_2_1>=0;y3_2_2>=0;y3_2_3>=0;y3_6_1>=0;y3_6_2>=0;y3_6_3>=0;
y4_6_1>=0;y4_6_2>=0;y4_6_3>=0;y4_6_4>=0;y4_5_1>=0;y4_5_2>=0;y4_5_3>=0;y4_5_4>=0;y4_1_1>=0;y4_1_2>=0;y4_1_3>=0;
y5_6_1>=0;y5_6_2>=0;y5_6_3>=0;y5_6_4>=0;y5_6_5>=0;y5_6_6>=0;y5_6_7>=0;y5_6_8>=0;y5_4_1>=0;y5_4_2>=0;y5_4_3>=0;y5_4_4>=0;y5_1_1>=0;y5_1_2>=0;y5_2>=0;y5_3>=0;
y1_6_1>=0;y1_6_2>=0;y1_6_3>=0;y1_6_4>=0;y1_4_1>=0;y1_4_2>=0;y1_4_3>=0;y1_5_1>=0;y1_5_2>=0;y1_2>=0;

y2_1<=1;y2_6_1<=1;y2_6_2<=1;y2_3_1<=1;y2_3_2<=1;y2_3_3<=1;y2_5<=1;
y3_5<=1;y3_2_1<=1;y3_2_2<=1;y3_2_3<=1;y3_6_1<=1;y3_6_2<=1;y3_6_3<=1;
y4_6_1<=1;y4_6_2<=1;y4_6_3<=1;y4_6_4<=1;y4_5_1<=1;y4_5_2<=1;y4_5_3<=1;y4_5_4<=1;y4_1_1<=1;y4_1_2<=1;y4_1_3<=1;
y5_6_1<=1;y5_6_2<=1;y5_6_3<=1;y5_6_4<=1;y5_6_5<=1;y5_6_6<=1;y5_6_7<=1;y5_6_8<=1;y5_4_1<=1;y5_4_2<=1;y5_4_3<=1;y5_4_4<=1;y5_1_1<=1;y5_1_2<=1;y5_2<=1;y5_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
% Identifying 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);


