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

# bus data

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};

# branch data

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 branch_b       {(l,k,m) in BRANCH} :=-branch_x[l,k,m]/(branch_r[l,k,m]^2+branch_x[l,k,m]^2);

# active power generation limits

param p_gen_upper;
param p_gen_lower;

# variables

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];

# auxiliar variables

var p_g {BUS}; # final active power generation, used to output data
var q_g {BUS}; # final reactive power generation, used to output data

var p_d {BRANCH}; # final active direct flow, used to output data
var q_d {BRANCH}; # final reactive direct flow, used to output data
var p_r {BRANCH}; # final active reverse flow, used to output data
var q_r {BRANCH}; # final reactive 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);

var G{(k,m) in YBUS} =
if(k == m) then (bus_g_shunt[k] + sum{(l,k,i) in BRANCH} branch_g[l,k,i]*branch_tap[l,k,i]^2
                                + sum{(l,i,k) in BRANCH} branch_g[l,i,k])
else if(k != m) then (sum{(l,k,m) in BRANCH} (-branch_g[l,k,m]*cos(branch_def[l,k,m])-branch_b[l,k,m]*sin(branch_def[l,k,m]))*branch_tap[l,k,m]
                     +sum{(l,m,k) in BRANCH} (-branch_g[l,m,k]*cos(branch_def[l,m,k])+branch_b[l,m,k]*sin(branch_def[l,m,k]))*branch_tap[l,m,k]);
 
var B{(k,m) in YBUS} =
if(k == m) then (bus_b_shunt[k] + sum{(l,k,i) in BRANCH} (branch_b[l,k,i]*branch_tap[l,k,i]^2 + branch_c[l,k,i]/2)
                                + sum{(l,i,k) in BRANCH} (branch_b[l,i,k]+branch_c[l,i,k]/2))
else if(k != m) then (sum{(l,k,m) in BRANCH} (branch_g[l,k,m]*sin(branch_def[l,k,m])-branch_b[l,k,m]*cos(branch_def[l,k,m]))*branch_tap[l,k,m]
                     +sum{(l,m,k) in BRANCH} (-branch_g[l,m,k]*sin(branch_def[l,m,k])-branch_b[l,m,k]*cos(branch_def[l,m,k]))*branch_tap[l,m,k]);

# important information

var reactive_generation =  sum {k in BUS} abs(bus_q_load[k] + sum{(k,m) in YBUS}
                          ((G[k,m]*sin(bus_angle[k]-bus_angle[m])-B[k,m]*cos(bus_angle[k]-bus_angle[m]))));

var inductive_generation =  sum {k in BUS} min(bus_q_load[k] + sum{(k,m) in YBUS} 
                          ((G[k,m]*sin(bus_angle[k]-bus_angle[m])-B[k,m]*cos(bus_angle[k]-bus_angle[m]))),0);

var capacitive_generation =  sum {k in BUS} max(bus_q_load[k] + sum{(k,m) in YBUS} 
                          ((G[k,m]*sin(bus_angle[k]-bus_angle[m])-B[k,m]*cos(bus_angle[k]-bus_angle[m]))),0);

# objectives

minimize active_power : sum{k in BUS : bus_type[k] == 2 || bus_type[k] == 3} 
   (bus_p_load[k] + sum{(k,m) in YBUS} 
    ((G[k,m]*cos(bus_angle[k]-bus_angle[m])+B[k,m]*sin(bus_angle[k]-bus_angle[m]))))^2;

minimize losses : sum{(l,k,m) in BRANCH} (branch_g[l,k,m]*( 2
-2*bus_voltage[k]*bus_voltage[m]*branch_tap[l,k,m]*cos(bus_angle[k]-bus_angle[m])));

minimize reactive_power : sum{k in BUS : bus_type[k] == 2 || bus_type[k] == 3} 
   (bus_q_load[k] + sum{(k,m) in YBUS} 
     ((G[k,m]*sin(bus_angle[k]-bus_angle[m])-B[k,m]*cos(bus_angle[k]-bus_angle[m]))))^2;

# constraints

subject to p_load {k in BUS : bus_type[k] == 0}:
   bus_p_gen[k] - bus_p_load[k] - sum{(k,m) in YBUS} 
                          ((G[k,m]*cos(bus_angle[k]-bus_angle[m])+B[k,m]*sin(bus_angle[k]-bus_angle[m]))) = 0;

subject to p_inj {k in BUS : bus_type[k] == 2 || bus_type[k] == 3}:
#   p_gen_lower*bus_p_gen[k] <= bus_p_load[k] + sum{(k,m) in YBUS} 
   0 <= bus_p_load[k] + sum{(k,m) in YBUS} 
               ((G[k,m]*cos(bus_angle[k]-bus_angle[m])+B[k,m]*sin(bus_angle[k]-bus_angle[m]))) <= p_gen_upper*bus_p_gen[k];


# data specifications

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  IEEE162a.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  IEEE162a.branch;

param p_gen_upper := 1.10;
param p_gen_lower := 0.90;

# data scaling and initialization

for{i in BUS} {
   let bus_voltage[i] := 1;
   let bus_angle[i] := 0;
   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;
  };


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 branch_tap[l,k,m] := 1;
  };

# fixed variables

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 superbasics_limit = 200 major_iterations = 100 timing=1";
option loqo_options "verbose=2 iterlim=200 timing=1";
option lancelot_options "timing=1";
option snopt_options "outlev = 1 superbasics_limit = 1000 timing=1"; 
solve active_power;

# calculates both active and reactive power generation

for{k in BUS} { 
  let p_g[k]  := bus_p_load[k] + sum{(k,m) in YBUS} 
                 ((G[k,m]*cos(bus_angle[k]-bus_angle[m])+B[k,m]*sin(bus_angle[k]-bus_angle[m])));


  let q_g[k]  := 0; 
}

# calculates both active and reactive direct and reverse fluxes

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

let p_d[l,k,m] := branch_g[l,k,m] 
-branch_g[l,k,m]*cos(bus_angle[k]-bus_angle[m]+branch_def[l,k,m])
-branch_b[l,k,m]*sin(bus_angle[k]-bus_angle[m]+branch_def[l,k,m]);

let q_d[l,k,m] := 0;

let p_r[l,k,m] := branch_g[l,k,m] 
-branch_g[l,k,m]*cos(bus_angle[k]-bus_angle[m]+branch_def[l,k,m])
+branch_b[l,k,m]*sin(bus_angle[k]-bus_angle[m]+branch_def[l,k,m]);


let q_r[l,k,m] := 0;
}

# generates output file


printf "Active Power Losses         %8.2f MW\n", losses*100 > aopf.stt;
printf "Total Active Generation     %8.2f MW\n", sum {k in BUS} p_g[k]*100 >> aopf.stt;
printf "Total Reactive Generation   %8.2f MVAR\n", reactive_generation*100 >> aopf.stt;
printf "Total Inductive Generation  %8.2f MVAR\n", inductive_generation*100 >> aopf.stt;
printf "Total Capacitive Generation %8.2f MVAR\n\n", capacitive_generation*100 >> aopf.stt;
printf "  #      Name    Voltage  Angle     PGen     QGen    PLoad    QLoad   Qshunt  To    PFlux    QFlux\n" >> aopf.stt;
printf "--------------------------------------------------------------------------------------------------\n" >> aopf.stt; 
for{i in BUS} {
printf "%4d %s %6.4f %6.2f %8.2f %8.2f %8.2f %8.2f %8.2f", i, bus_name[i], bus_voltage[i], bus_angle[i]*180/3.14159,
p_g[i]*100, q_g[i]*100, bus_p_load[i]*100, bus_q_load[i]*100, bus_voltage[i]^2*bus_b_shunt[i]*100 >> aopf.stt;

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

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

for{(l,k,i) in BRANCH} {
printf "%75s %4d %8.2f %8.2f\n", "", k, p_r[l,k,i]*100, q_r[l,k,i]*100 >> aopf.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], 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],
1/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;