set BUS;    # set of buses
set BRANCH within {1..4000} cross BUS cross BUS; # set of branches

param bus_type       {BUS};
param bus_name       {BUS} symbolic;
param bus_voltage    {BUS};
param bus_angle0     {BUS};
param bus_p_gen      {BUS};
param bus_q_gen      {BUS};
param bus_q_min      {BUS};
param bus_q_max      {BUS};
param bus_p_load     {BUS};
param bus_q_load     {BUS};
param bus_g_shunt    {BUS};
param bus_b_shunt    {BUS};
param bus_b_shunt_min{BUS};
param bus_b_shunt_max{BUS};
param bus_b_dispatch {BUS};
param bus_area       {BUS};

param branch_from   {BRANCH};
param branch_to     {BRANCH};
param branch_type   {BRANCH};
param branch_r      {BRANCH};
param branch_x      {BRANCH};
param branch_c      {BRANCH};
param branch_tap    {BRANCH};
param branch_tap_min{BRANCH};
param branch_tap_max{BRANCH};
param branch_def0   {BRANCH};
param branch_def_min{BRANCH};
param branch_def_max{BRANCH};


param branch_g {(l,k,m) in BRANCH} := branch_r[l,k,m]/(branch_r[l,k,m]^2+branch_x[l,k,m]^2);
param losses;
param dist_losses{BUS};
param p_gen_upper;

var bus_angle   {i in BUS};
var branch_def  {(l,k,m) in BRANCH} >= branch_def_min[l,k,m], <= branch_def_max[l,k,m];

var p_d {BRANCH}; # final active direct flow, used to output data
var p_r {BRANCH}; # final active reverse flow, used to output data

# matrix YBUS

set YBUS := setof{i in BUS} (i,i) union 
setof {(l,k,m) in BRANCH} (k,m) union
setof {(l,k,m) in BRANCH} (m,k);

param B{(k,m) in YBUS} :=
if(k == m) then ( sum{(l,k,i) in BRANCH} (1/branch_x[l,k,i]) 
                + sum{(l,i,k) in BRANCH} (1/branch_x[l,i,k]))
else if(k != m) then (sum{(l,k,m) in BRANCH} 1/branch_x[l,k,m]
                     +sum{(l,m,k) in BRANCH} 1/branch_x[l,m,k]);


# calculates both active and reactive power generation

var p_g {k in BUS} = bus_p_load[k] + dist_losses[k] + sum{(k,m) in YBUS} 
                 (B[k,m]*(bus_angle[k]-bus_angle[m]+sum{(l,k,m) in BRANCH} branch_def[l,k,m]+sum{(l,m,k) in BRANCH} -branch_def[l,m,k]));
 

minimize active_power : sum {k in BUS : bus_type[k] == 2 || bus_type[k] == 3}
(bus_p_load[k]+sum{(k,m) in YBUS} B[k,m]*(bus_angle[k]-bus_angle[m]+sum{(l,k,m) in BRANCH} branch_def[l,k,m]+sum{(l,m,k) in BRANCH} -branch_def[l,m,k]));
 
subject to p_load {k in BUS : bus_type[k] == 0}:
   bus_p_gen[k] - bus_p_load[k] - dist_losses[k] - sum{(k,m) in YBUS} 
                        (B[k,m]*(bus_angle[k]-bus_angle[m]+sum{(l,k,m) in BRANCH} branch_def[l,k,m]+sum{(l,m,k) in BRANCH} -branch_def[l,m,k] )) = 0;

subject to p_gen {k in BUS : bus_type[k] == 2 || bus_type[k] == 3}:
  0  <= bus_p_load[k] + dist_losses[k] + sum{(k,m) in YBUS} (B[k,m]*(bus_angle[k]-bus_angle[m]+
sum{(l,k,m) in BRANCH} branch_def[l,k,m]+sum{(l,m,k) in BRANCH} -branch_def[l,m,k] )) <= p_gen_upper*bus_p_gen[k];

data;

param: BUS: bus_type bus_name bus_voltage bus_angle0 bus_p_gen bus_q_gen
            bus_q_min bus_q_max bus_p_load bus_q_load bus_g_shunt bus_b_shunt
            bus_b_shunt_min bus_b_shunt_max bus_b_dispatch bus_area := 
include IEEE662a.bus;

param: BRANCH: branch_type branch_r branch_x branch_c
               branch_tap branch_tap_min branch_tap_max branch_def0 
               branch_def_min branch_def_max :=
include IEEE662a.branch;


for{i in BUS} {
   let bus_angle[i] := bus_angle0[i]*3.14159/180;
   let bus_p_gen[i] := bus_p_gen[i]/100;
   let bus_q_gen[i] := bus_q_gen[i]/100;
   let bus_q_min[i] := bus_q_min[i]/100;
   let bus_q_max[i] := bus_q_max[i]/100;
   let bus_p_load[i] := bus_p_load[i]/100;
   let bus_q_load[i] := bus_q_load[i]/100;
   let dist_losses[i] :=0;
  };

for{(l,k,m) in BRANCH} {
   let branch_def[l,k,m] := -branch_def0[l,k,m]*3.14159/180;
   let branch_def_min[l,k,m] := branch_def_min[l,k,m]*3.14159/180;
   let branch_def_max[l,k,m] := branch_def_max[l,k,m]*3.14159/180;  
};

let p_gen_upper := 1.10;

fix {i in BUS : bus_type[i] == 3} bus_angle[i]; # slack angle fixed
fix {(l,k,m) in BRANCH : branch_type[l,k,m] !=4} branch_def[l,k,m]; # no phase shifters on branch_type different of 4

option minos_options "summary_file=6 timing=1";
option loqo_options "verbose=2 timing=1";
option lancelot_options "timing=1";
option snopt_options "timing=1";
solve;

let losses := sum {(l,k,m) in BRANCH} branch_g[l,k,m]*(bus_angle[k]-bus_angle[m]+branch_def[l,k,m])^2;

param active_load := sum {i in BUS} abs(bus_p_load[i]); 

for{i in BUS} {
let dist_losses[i] := abs(bus_p_load[i])*losses/active_load;
};

solve;

# calculates both active and reactive direct and reverse fluxes

for{(l,k,m) in BRANCH} {

let p_d[l,k,m] := (bus_angle[k]-bus_angle[m]+branch_def[l,k,m])/branch_x[l,k,m];

let p_r[l,k,m] := - p_d[l,k,m];
}

# generates output file

printf "Losses            %8.2f MW\n", losses*100 > dcopf.stt;
printf "Active Generation %8.2f MW\n", sum {k in BUS} p_g[k]*100 >> dcopf.stt;
printf "  #      Name    Voltage  Angle     PGen    PLoad  To    PFlux\n" >> dcopf.stt;
printf "--------------------------------------------------------------\n" >> dcopf.stt; 
for{i in BUS} {
printf "%4d %s %6.4f %6.2f %8.2f %8.2f", i, bus_name[i], bus_voltage[i], bus_angle[i]*180/3.14159,
p_g[i]*100, bus_p_load[i]*100 >> dcopf.stt;

printf " ------------\n" >> dcopf.stt;

for{(l,i,m) in BRANCH} 
printf "%48s %4d %8.2f\n", "", m, p_d[l,i,m]*100 >> dcopf.stt;

for{(l,k,i) in BRANCH } 
printf "%48s %4d %8.2f\n", "", k, p_r[l,k,i]*100 >> dcopf.stt;
}


# generates output for next input

for{i in BUS}
printf "%4d %1d %c%s%c %8.4f %8.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %8.4f %8.4f %8.4f %8.4f %1d %2d\n", 
i, bus_type[i], '"', bus_name[i], '"',bus_voltage[i], 0, #bus_angle[i]*180/3.14159,
p_g[i]*100, bus_q_gen[i]*100, bus_q_min[i]*100, bus_q_max[i]*100, bus_p_load[i]*100, 
bus_q_load[i]*100, bus_g_shunt[i], bus_b_shunt[i], bus_b_shunt_min[i],
bus_b_shunt_max[i], bus_b_dispatch[i], bus_area[i] > out.bus;

for{(l,k,m) in BRANCH}
printf "%4d %4d %4d %1d %10.6f %10.6f %10.6f %6.4f %6.4f %6.4f %6.2f %6.2f %6.2f\n",
l, k, m, branch_type[l,k,m], branch_r[l,k,m], branch_x[l,k,m], branch_c[l,k,m],
branch_tap[l,k,m], branch_tap_min[l,k,m], branch_tap_max[l,k,m], 
-branch_def[l,k,m]*180/3.14159, branch_def_min[l,k,m]*180/3.14159,
branch_def_max[l,k,m]*180/3.14159 > out.branch;