########
# Sets #
########

set BUS;    # set of buses
set BRANCH within {1..4000} cross BUS cross BUS; # set of branches
set GEN;    # set of GENs
set LOAD    ordered;# set of loads
set PHASES; # set of phases {a, b, c}
set LINES  ordered;  # set of lines
set ENDS;   # set of ends {s, r} s: sending end, r: receiving end
set LAYERS; # set of layers {x, y} x: real part, y: imaginary part
set GENS;    # set of generators
set XSL;
set XSLVAR default {}; # set of variable XSL
set OLTC default {};

param Pfactor ;
param Qfactor ;
param Gfactor ;
##############
# Parameters #
##############

# general constants
param pi:=3.141592654; 
# node data
param node_type{n in BUS, p in PHASES}; # A node is defined by its bus 	
	# number and phase number.

param LtoLE{PHASES, PHASES}; # matrix to conver line-to-ground voltages
							# to line-to-line values.
param LtoLI{PHASES, PHASES}; # matrix to conver line-to-ground currents
							# to line-to-line values.


# line data
#~~~~~~~~~~~~
# this can be used to model lines, transformer, LTC, and switch.
# For every GEN there is a line representing X_SL
# (don't forget to adjust the MVAbase!)


# voltage at each node. a node is defined by its bus number n and phase p.
#~~~~~~~~~~~~~~~~~~~~~~~
var Ex {n in BUS, p in PHASES} >= -2.0 <= 2.0; 
# real voltage at bus n and phase p
var Ey {n in BUS, p in PHASES} >= -2.0 <= 2.0; 
# imaginary Voltage at bus n and phase p
var Emag{n in BUS, p in PHASES} 
    = sqrt(Ex[n,p]^2 + Ey[n,p]^2); # magnitude of voltages
var Eangle{n in BUS, p in PHASES} 
    = atan2(Ey[n, p], Ex[n, p])*180/pi; # angle of voltages
	
var ExLL {n in BUS, p in PHASES} = sum {pa in PHASES} 
                                   ( LtoLE[p, pa] *Ex[n,pa])/sqrt(3);
var EyLL {n in BUS, p in PHASES} = sum {pa in PHASES} 
                                   ( LtoLE[p, pa] *Ey[n,pa])/sqrt(3);
var ELLmag{n in BUS, p in PHASES} 
    = sqrt(ExLL[n,p]^2 + EyLL[n,p]^2); # magnitude of voltages

var ELLangle{n in BUS, p in PHASES} = atan2(EyLL[n, p], 
									        ExLL[n, p])*180/pi; 
								# angle of voltages


## Constant impedance load
# current equation
#~~~~~~~~~~~~~~~~~~~

# imag gen's voltage
                                        # connected at bus n phase p

										# generator data

# var gen_p{g in GENS, n in BUS, p in PHASES}; 
# real generation at bus n
	# # phase p
# var gen_q{g in GENS, n in BUS, p in PHASES}; # imag generation at bus n
	# # phase p

	
# current balance equation
#~~~~~~~~~~~~~~~~~~~~~~~~~~~
var I_gx {n in BUS, p in PHASES}; # sum of real current of ALL generators
        # connected at bus n phase p
var I_gy {n in BUS, p in PHASES}; # sum imag current of ALL generators

# General quantities of interest
# Real and Imaginary power consumed at substation
# var Psub = sum {p in PHASES} line_ps['2_3', p]; # the first line is '2_3'; '1_2' is OLTC
# var Qsub = sum {p in PHASES} line_qs['2_3', p]; # the first line is '2_3'; '1_2' is OLTC

################################
# LINE Model
################################
# include LINE.mod
# This is the line model

param line_name{LINES} symbolic;
param line_type{LINES}; # types:
                # 1. Lines
                # 2. transformers
                # 3. Load tap changers
                # 4. gen lines
param line_number{LINES}; # number of a line
param line_Ax{LINES, PHASES, PHASES}; # real part of A matrix
param line_Ay{LINES, PHASES, PHASES}; # imaginary part of A matrix
param line_Bx{LINES, PHASES, PHASES}; # real part of B matrix
param line_By{LINES, PHASES, PHASES}; # imaginary part of B matrix
param line_C{LINES, PHASES, PHASES}; # real part of C matrix
param line_Bmag {li in LINES, p1 in PHASES, p2 in PHASES} 
      := sqrt ( line_Bx[li, p1, p2]^2 + line_By[li, p1, p2]^2) ;
param line_R {li in LINES, p1 in PHASES, p2 in PHASES
             : line_Bmag[li, p1, p2] != 0} 
			 := line_Bx[li, p1, p2]/line_Bmag[li, p1, p2] ;
##param line_Cx{LINES, PHASES, PHASES}; # real part of C matrix
##param line_Cy{LINES, PHASES, PHASES}; # imag part of C matrix
param line_D{LINES, PHASES, PHASES}; # real part of D matrix
##param line_Dx{LINES, PHASES, PHASES}; # real part of D matrix
##param line_Dy{LINES, PHASES, PHASES}; # imag part of D matrix
param line_sbus{li in LINES}; # a line's sending bus
param line_rbus{li in LINES}; # a line's receiving bus


# variable for series component
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
var I_lxr { li in LINES, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a line's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lxs { li in LINES, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a line's
        # e end connected to bus n phase p. e in {receiving, sending}

var I_lyr { li in LINES, p in PHASES} >= -2.0 <= 2.0;
		# Imag current of a line's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lys { li in LINES, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a line's
        # e end connected to bus n phase p. e in {receiving, sending}

var I_lr2 { li in LINES, p in PHASES} =  I_lxr [li, p]^2 + I_lyr [li, p]^2 ;
var I_lrmag { li in LINES, p in PHASES} 
	=  sqrt (I_lxr [li, p]^2 + I_lyr [li, p]^2 + 1e-48);
var I_lrangle { li in LINES, p in PHASES} 
	=  atan2(I_lyr [li, p], I_lxr [li, p])*180/pi ;

# var line_loss { li in LINES, p in PHASES} = 
	# sum {pa in PHASES : line_Bmag[li, p, pa] != 0} 
	# ( line_Bx[li, p, pa] * I_lr2[li, pa]/line_Bmag[li, p, pa] ) ;

var line_loss { li in LINES, p in PHASES} = 
	sum {pa in PHASES : line_Bmag[li, p, pa] != 0 and p=pa} 
	( line_Bx[li, p, pa] * I_lr2[li, pa]) ;

var line_pr { li in LINES, p in PHASES} =
     Ex[line_rbus[li],p]*I_lxr[li, p] 
	+ Ey[line_rbus[li],p]*I_lyr[li, p] ; 

var line_ps { li in LINES, p in PHASES} =
     Ex[line_sbus[li],p]*I_lxs[li, p] 
	+ Ey[line_sbus[li],p]*I_lys[li, p] ; 

var line_qs { li in LINES, p in PHASES} =
    - Ex[line_sbus[li],p]*I_lys[li, p] 
	+ Ey[line_sbus[li],p]*I_lxs[li, p] ; 

var line_qr { li in LINES, p in PHASES} =
    - Ex[line_rbus[li],p]*I_lyr[li, p] 
	+ Ey[line_rbus[li],p]*I_lxr[li, p] ; 
	
var line_loss2 { li in LINES, p in PHASES} = 
	# abs(line_ps[li,p, ti]-line_pr[li,p, ti]); abs is discontinous function at origin. I think this is why it is difficult to solve the line loss obj.
	line_Bx[li,p,p]*I_lrmag[li,p]^2 ;

# var line_loss2q { li in LINES, p in PHASES} = 
	# abs(line_qs[li,p]-line_qr[li,p]);

	
var all_line_loss = sum { li in LINES, p in PHASES}
                     line_loss2[li, p];

# var all_line_lossq = sum { li in LINES, p in PHASES}
                     # line_loss2q[li, p];
					 
var line_total_loss = sum { li in LINES, p in PHASES} line_loss[li, p];
param line_xsl{li in LINES : line_type[li]=4} >=0;

# General quantities of interest
# # Real and Imaginary power consumed at substation
var Psub = sum {p in PHASES} line_ps['149_1', p]; # the first line is '2_3'; '1_2' is OLTC
var Qsub = sum {p in PHASES} line_qs['149_1', p]; # the first line is '2_3'; '1_2' is OLTC


################################
# End of LINE Model
################################
################################
# LOAD Model
################################
#~~~~~~~~~~~~
param load_type{LOAD}; # type can be constant power, impedance or currrent.
# Load types are:
# 1. Constant power
# 2. Constant impedance
# 3. Constant current
param load_bus{LOAD}; # The node to which the load is connected.
param load_phase{LOAD}; # The phase to which the load is connected.
param load_p{L in LOAD, n in BUS, p in PHASES}; # the real power constant
        # load. ===> for constant power load only
param load_q{L in LOAD, n in BUS, p in PHASES}; # the reactive power
        # constant load. ===> for constant power load only
param load_pD{L in LOAD, n in BUS, p in PHASES}; # the real power constant
        # load. ===> for constant power load only
param load_qD{L in LOAD, n in BUS, p in PHASES}; # the reactive power
param load_r{L in LOAD, n in BUS, p in PHASES}; # the resistance of
        # a constant impedance load.
param load_x{L in LOAD, n in BUS, p in PHASES}; # the inductance of
        # a constant impedance load.
param load_rD{L in LOAD, n in BUS, p in PHASES}; # the resistance of
        # a constant impedance load.
param load_xD{L in LOAD, n in BUS, p in PHASES}; # the inductance of
        # a constant impedance load.
param load_Imag{L in LOAD, n in BUS, p in PHASES};  # the current
        #magnitude of a constant current load
param load_PFangle {L in LOAD, n in BUS, p in PHASES}; # the power
		# factor angle of a constant current load
param load_ImagD {L in LOAD, n in BUS, p in PHASES};  # the current
        #magnitude of a constant current load
param load_PFangleD {L in LOAD, n in BUS, p in PHASES}; # the power
		# factor angle of a constant current load
param load_connection{LOAD};


# Load connections are:
# 1 wye
# 2 delta

# variable for load
#~~~~~~~~~~~~~~~~~~~~

## Constant power load
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
var E_Lx {L in LOAD, p in PHASES: load_type[L]=1} = Ex[load_bus[L],p]; 
# Real voltage of a load
        # connected to bus n phase p
var E_Ly {L in LOAD, p in PHASES: load_type[L]=1} = Ex[load_bus[L],p]; 
# Imaginary voltage of a load
        # connected to bus n phase p
var I_Lx {L in LOAD, p in PHASES: load_type[L]=1  } = (  
              load_p[L, load_bus[L], p]*Pfactor * Ex[load_bus[L], p] 
		    + load_q[L, load_bus[L], p]*Qfactor * Ey[load_bus[L], p]  ) 
			  / Emag[load_bus[L], p]; # Real current of a load
        # connected to bus n phase p
var I_Ly {L in LOAD, p in PHASES: load_type[L]=1 } = (  
              load_p[L, load_bus[L], p]*Pfactor * Ey[load_bus[L], p] 
            - load_q[L, load_bus[L], p]*Qfactor * Ex[load_bus[L], p]  ) 
			/ Emag[load_bus[L], p]; # Imag voltage of a load
        # connected to bus n phase p
		
var I_LxLL {L in LOAD, p in PHASES: load_type[L]=1  } = (  
			  load_pD[L, load_bus[L], p]*Pfactor * ExLL[load_bus[L], p] 
		    + load_qD[L, load_bus[L], p]*Qfactor * EyLL[load_bus[L], p]  ) 
			  / (ELLmag[load_bus[L], p]);#*sqrt(3));; # Real current of a load
        # connected to bus n phase p

var I_LyLL {L in LOAD, p in PHASES: load_type[L]=1 } = (  
              load_pD[L, load_bus[L], p]*Pfactor * EyLL[load_bus[L], p] 
            - load_qD[L, load_bus[L], p]*Qfactor * ExLL[load_bus[L], p]  ) 
			/ (ELLmag[load_bus[L], p]);#*sqrt(3)); # Imag voltage of a load
        # connected to bus n phase p

# for a delta connected CS load: corresponding real phase current
var I_LxD {L in LOAD, p in PHASES: load_type[L]=1 } 
		= sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL[L,pa]/sqrt(3));
		# =  sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL[L,pa]);
# for a delta connected CS load: corresponding real phase current
var I_LyD {L in LOAD, p in PHASES: load_type[L]=1} 
         = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL[L,pa]/sqrt(3));
		 # = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL[L,pa]);

# wye-connected CS LOAD		 
var I_Lmag {L in LOAD, p in PHASES: load_type[L]=1} 
			= sqrt(I_Lx[L,p]^2 + I_Ly[L,p]^2);
var I_Langle {L in LOAD, p in PHASES: load_type[L]=1 } 
			= atan2(I_Ly[L,p], I_Lx[L,p])*180/pi;
var I_LLLmag {L in LOAD, p in PHASES : load_type[L]=1} 
			= sqrt(I_LxLL[L,p]^2 + I_LyLL[L,p]^2);
var I_LLLangle {L in LOAD, p in PHASES : load_type[L]=1} 
			= atan2(I_LyLL[L,p], I_LxLL[L,p])*180/pi;
var I_LDmag {L in LOAD, p in PHASES : load_type[L]=1} 
			= sqrt(I_LxD[L,p]^2 + I_LyD[L,p]^2);
var I_LDangle {L in LOAD, p in PHASES: load_type[L]=1 } 
			= atan2(I_LyD[L,p], I_LxD[L,p])*180/pi;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~			
# Constant Impedance Load
# Real current
var I_Lx_Imp {L in LOAD, p in PHASES: load_type[L]=2}
         = ( Ex[load_bus[L], p]*load_r[L , load_bus[L], p]
		     + Ey[load_bus[L], p] *load_x[L , load_bus[L], p] )
		   /( load_r[L , load_bus[L], p]^2 + load_x[L , load_bus[L], p]^2 ) ;
		   
		   
#Imaginary current	
var I_Ly_Imp {L in LOAD, p in PHASES: load_type[L]=2}
         = ( - Ex[load_bus[L], p] * load_x[L , load_bus[L], p]
		     + Ey[load_bus[L], p] * load_r[L , load_bus[L], p] )
		   /( load_r[L , load_bus[L], p]^2 + load_x[L , load_bus[L], p]^2 ) ;

# line-to-line current
var I_LxLL_Imp {L in LOAD, p in PHASES: load_type[L]=2  } 
         = ( ExLL[load_bus[L], p]*load_rD[L , load_bus[L], p]
		     + EyLL[load_bus[L], p] *load_xD[L , load_bus[L], p] ) 
		   /( load_rD[L , load_bus[L], p]^2 + load_xD[L , load_bus[L], p]^2 ) ;
		   	
		 

var I_LyLL_Imp {L in LOAD, p in PHASES: load_type[L]=2 }
         = ( - ExLL[load_bus[L], p] * load_xD[L , load_bus[L], p]
		     + EyLL[load_bus[L], p] * load_rD[L , load_bus[L], p] )
		   /( load_rD[L , load_bus[L], p]^2 + load_xD[L , load_bus[L], p]^2 ) ;
		   
		   
# for a delta connected CI load: corresponding real phase current
var I_LxD_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
		= sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_Imp[L,pa]/sqrt(3));
		# = - sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_Imp[L,pa]);
# adding a negative sign lead to potive drop voltage over line 78. 
# But why is that so? please carefully examine this!
# for a delta connected CS load: corresponding real phase current
var I_LyD_Imp {L in LOAD, p in PHASES: load_type[L]=2} 
         = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_Imp[L,pa]/sqrt(3));
		 # = - sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_Imp[L,pa]);
		 
		 
# wye-connected CS LOAD		 
var I_Lmag_Imp {L in LOAD, p in PHASES : load_type[L]=2 } 
			= sqrt(I_Lx_Imp[L,p]^2 + I_Ly_Imp[L,p]^2);
var I_Langle_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
			= atan2(I_Ly_Imp[L,p], I_Lx_Imp[L,p])*180/pi;
var I_LLLmag_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
			= sqrt(I_LxLL_Imp[L,p]^2 + I_LyLL_Imp[L,p]^2);
var I_LLLangle_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
			= atan2(I_LyLL_Imp[L,p], I_LxLL_Imp[L,p])*180/pi;
var I_LDmag_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
			= sqrt(I_LxD_Imp[L,p]^2 + I_LyD_Imp[L,p]^2);
var I_LDangle_Imp {L in LOAD, p in PHASES: load_type[L]=2 } 
			= atan2(I_LyD_Imp[L,p], I_LxD_Imp[L,p])*180/pi;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Constant Current Load
# Real current
var I_Lx_Cur {L in LOAD, p in PHASES: load_type[L]=3}
         = load_Imag[L, load_bus[L], p]*cos( atan2( Ey[load_bus[L], p], Ex[load_bus[L], p] ) - load_PFangle[L, load_bus[L], p] ) ;
		   
		   
#Imaginary current	
var I_Ly_Cur {L in LOAD, p in PHASES: load_type[L]=3}
         = load_Imag[L, load_bus[L], p]*sin( atan2( Ey[load_bus[L], p], Ex[load_bus[L], p] ) - load_PFangle[L, load_bus[L], p] ) ;

# line-to-line current
var I_LxLL_Cur {L in LOAD, p in PHASES: load_type[L]=3  } 
         = load_ImagD[L, load_bus[L], p]*cos( atan2( EyLL[load_bus[L], p], ExLL[load_bus[L], p] ) - load_PFangleD[L, load_bus[L], p] );#/sqrt(3);
# Why should it be divided by sqrt(3)? Because its delta connected.
		   	
		 

var I_LyLL_Cur {L in LOAD, p in PHASES: load_type[L]=3 }
         = load_ImagD[L, load_bus[L], p]*sin( atan2( EyLL[load_bus[L], p], ExLL[load_bus[L], p] ) - load_PFangleD[L, load_bus[L], p] );#/sqrt(3); 
		   
		   
# for a delta connected CI load: corresponding real phase current
var I_LxD_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
		= sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_Cur[L,pa]/sqrt(3));
		# =  sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_Cur[L,pa]);
# for a delta connected CS load: corresponding real phase current
var I_LyD_Cur {L in LOAD, p in PHASES: load_type[L]=3} 
         = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_Cur[L,pa]/sqrt(3));
		 # =  sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_Cur[L,pa]);
		 
		 
# wye-connected CI LOAD		 
var I_Lmag_Cur {L in LOAD, p in PHASES : load_type[L]=3 } 
			= sqrt(I_Lx_Cur[L,p]^2 + I_Ly_Cur[L,p]^2);
var I_Langle_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
			= atan2(I_Ly_Cur[L,p], I_Lx_Cur[L,p])*180/pi;
var I_LLLmag_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
			= sqrt(I_LxLL_Cur[L,p]^2 + I_LyLL_Cur[L,p]^2);
var I_LLLangle_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
			= atan2(I_LyLL_Cur[L,p], I_LxLL_Cur[L,p])*180/pi;
var I_LDmag_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
			= sqrt(I_LxD_Cur[L,p]^2 + I_LyD_Cur[L,p]^2);
var I_LDangle_Cur {L in LOAD, p in PHASES: load_type[L]=3 } 
			= atan2(I_LyD_Cur[L,p], I_LxD_Cur[L,p])*180/pi;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## DG as negative constant power load
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # connected to bus n phase p
var I_Lx_dg {L in LOAD, p in PHASES: load_type[L]=4  } = (  
              load_p[L, load_bus[L], p]*Gfactor* Ex[load_bus[L], p] 
		    + load_q[L, load_bus[L], p]*Gfactor* Ey[load_bus[L], p]  ) 
			  / Emag[load_bus[L], p]; # Real current of a load
        # connected to bus n phase p
		
var I_Ly_dg {L in LOAD, p in PHASES: load_type[L]=4 } = (  
              load_p[L, load_bus[L], p]*Gfactor* Ey[load_bus[L], p] 
            - load_q[L, load_bus[L], p]*Gfactor* Ex[load_bus[L], p]  ) 
			/ Emag[load_bus[L], p]; # Imag voltage of a load
        # connected to bus n phase p
		
var I_LxLL_dg {L in LOAD, p in PHASES: load_type[L]=4  } = (  
			  load_pD[L, load_bus[L], p]*Gfactor* ExLL[load_bus[L], p] 
		    + load_qD[L, load_bus[L], p]*Gfactor* EyLL[load_bus[L], p]  ) 
			  / (ELLmag[load_bus[L], p]);#*sqrt(3));; # Real current of a load
        # connected to bus n phase p

var I_LyLL_dg {L in LOAD, p in PHASES: load_type[L]=4 } = (  
              load_pD[L, load_bus[L], p]*Gfactor* EyLL[load_bus[L], p] 
            - load_qD[L, load_bus[L], p]*Gfactor* ExLL[load_bus[L], p]  ) 
			/ (ELLmag[load_bus[L], p]);#*sqrt(3)); # Imag voltage of a load
        # connected to bus n phase p

# for a delta connected CS load: corresponding real phase current
var I_LxD_dg {L in LOAD, p in PHASES: load_type[L]=4 } 
		= sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_dg[L,pa]/sqrt(3));
		# =  sum {pa in PHASES} ( LtoLI[p, pa] *I_LxLL_dg[L,pa]);
# for a delta connected CS load: corresponding real phase current
var I_LyD_dg {L in LOAD, p in PHASES: load_type[L]=4} 
         = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_dg[L,pa]/sqrt(3));
		 # = sum {pa in PHASES}( LtoLI[p, pa]*I_LyLL_dg[L,pa]);

# wye-connected CS LOAD		 
var I_Lmag_dg {L in LOAD, p in PHASES: load_type[L]=4} 
			= sqrt(I_Lx_dg[L,p]^2 + I_Ly_dg[L,p]^2);
var I_Langle_dg {L in LOAD, p in PHASES: load_type[L]=4 } 
			= atan2(I_Ly_dg[L,p], I_Lx_dg[L,p])*180/pi;
var I_LLLmag_dg {L in LOAD, p in PHASES : load_type[L]=4} 
			= sqrt(I_LxLL_dg[L,p]^2 + I_LyLL_dg[L,p]^2);
var I_LLLangle_dg {L in LOAD, p in PHASES : load_type[L]=4} 
			= atan2(I_LyLL_dg[L,p], I_LxLL_dg[L,p])*180/pi;
var I_LDmag_dg {L in LOAD, p in PHASES : load_type[L]=4} 
			= sqrt(I_LxD_dg[L,p]^2 + I_LyD_dg[L,p]^2);
var I_LDangle_dg {L in LOAD, p in PHASES: load_type[L]=4 } 
			= atan2(I_LyD_dg[L,p], I_LxD_dg[L,p])*180/pi;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

################################
# End of LOAD Model
################################
################################
# GEN Model
################################
# include GEN.mod;
# model for generator


# generator data
# param gen_bus{GENS}; # the bus at which generator is connected.
param gen_P{GENS}; #
param gen_Q{GENS}; #
param gen_type{GENS}; # Type of generator.

# GEN Generator data
#~~~~~~~~~~~~~~~~~~~~~
param gen_line{GEN}; # GEN's line
param gen_gen{GEN}; # GEN's generator
param gen_qcrated{GEN}; # GEN's capacitive reactive rated power
param gen_qlrated{GEN}; # GEN's inductive real rated power
#param gen_vref{GEN}; # GEN's reference voltage
param gen_epsilon{GEN}; # GEN's epsilon
param gen_bus{GEN}; # the bus to which an GEN is connected.
param gen_bus_cont{GEN}; # the bus which is controlled by GEN



# variable for GEN
#~~~~~~~~~~~~~~~~~~~
# var gen_Vref{g in GENS};
var gen_vref{s in GEN} >= 0.95 <=1.05;# GEN's reference voltage
# var gen_XSL{g in GENS};
 # param gen_B {s in GEN};
#= gen_B[s]*Emag[gen_bus[s], p];
var gen_Ix {s in GEN, p in PHASES} >= -2.0 <= 2.0;
# The generator real current
				 
var gen_Iy {s in GEN, p in PHASES} >= -2.0 <= 2.0;
# The generator imag current
				 
# the generator active power
var gen_p {s in GEN, p in PHASES} = 
	  Ex[gen_bus[s],p]*gen_Ix[s, p] 
	+ Ey[gen_bus[s],p]*gen_Iy[s, p]; 
# the generator reactive power
var gen_q {s in GEN, p in PHASES} =
	- Ex[gen_bus[s],p]*gen_Iy[s, p] 
	+ Ey[gen_bus[s],p]*gen_Ix[s, p];	
# the real component of positive sequence voltage
var   Ex1{s in GEN} = (   Ex[gen_bus_cont[s],1] 
					  + Ex[gen_bus_cont[s],2]*cos(2*pi/3)
                      - Ey[gen_bus_cont[s],2]*sin(2*pi/3) 
					  + Ex[gen_bus_cont[s],3]*cos(4*pi/3)
					  - Ey[gen_bus_cont[s],3]*sin(4*pi/3) 
					 )/3;
# the imag component of positive sequence voltage
var Ey1{s in GEN} = (   Ey[gen_bus_cont[s],1] 
					  + Ey[gen_bus_cont[s],2]*cos(2*pi/3)
                      + Ex[gen_bus_cont[s],2]*sin(2*pi/3) 
					  + Ey[gen_bus_cont[s],3]*cos(4*pi/3)
					  + Ex[gen_bus_cont[s],3]*sin(4*pi/3) 
					 )/3	;
# the magnitude positive sequence voltage					 
var E1mag {s in GEN} = sqrt( Ex1[s]^2 + Ey1[s]^2);

# Variable for generators
#~~~~~~~~~~~~~~~~~~~~~~~~~~
var E_gx {g in GEN, p in PHASES} = Ex[gen_bus[g],p]; 
# real gen's voltage \# connected at bus n phase p

var E_gy {g in GEN, p in PHASES} = Ey[gen_bus[g],p]; 
################################
# End of GEN Model
################################
################################
# XSLVAR Model requires gen model
################################
# include XSLVAR.mod;
# This is the xslvar model
set XSLVAR_TAPS={0.0, 0.01, 0.02, 0.03, 0.04, 0.05};
set XSLVAR_TAPS_BIN={0, 1};

param xslvar_By_lower_bound {XSLVAR, PHASES};
param xslvar_name{XSLVAR} symbolic;
param xslvar_type{XSLVAR}; # types:
                # 1. Lines
                # 2. transformers
                # 3. Load tap changers
                # 4. gen xslvars
param xslvar_number{XSLVAR} >=0; # number of a xslvar
param xslvar_Ax{XSLVAR, PHASES, PHASES}; # real part of A matrix
param xslvar_Ay{XSLVAR, PHASES, PHASES}; # imaginary part of A matrix
param xslvar_Bx{XSLVAR, PHASES, PHASES}; # real part of B matrix
var xslvar_By_bin{XSLVAR, PHASES} in XSLVAR_TAPS_BIN;
var xslvar_By{XSLVAR, PHASES, PHASES} >=0 <=0.05;


var xslvar_xsl{XSLVAR} default 0;
 # imaginary part of B matrix
param xslvar_C{XSLVAR, PHASES, PHASES}; # real part of C matrix
# param xslvar_Bmag {li in XSLVAR, p1 in PHASES, p2 in PHASES} 
      # = sqrt ( xslvar_Bx[li, p1, p2]^2 + xslvar_By[li, p1, p2]^2) ;
# param xslvar_R {li in XSLVAR, p1 in PHASES, p2 in PHASES
             # : xslvar_Bmag[li, p1, p2] != 0} 
			 # := xslvar_Bx[li, p1, p2]/xslvar_Bmag[li, p1, p2] ;
##param xslvar_Cx{XSLVAR, PHASES, PHASES}; # real part of C matrix
##param xslvar_Cy{XSLVAR, PHASES, PHASES}; # imag part of C matrix
param xslvar_D{XSLVAR, PHASES, PHASES}; # real part of D matrix
##param xslvar_Dx{XSLVAR, PHASES, PHASES}; # real part of D matrix
##param xslvar_Dy{XSLVAR, PHASES, PHASES}; # imag part of D matrix
param xslvar_sbus{li in XSLVAR}; # a xslvar's sending bus
param xslvar_rbus{li in XSLVAR}; # a xslvar's receiving bus

# variable for series component
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
var I_lxr_xslvar { li in XSLVAR, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a xslvar's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lxs_xslvar { li in XSLVAR, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a xslvar's
        # e end connected to bus n phase p. e in {receiving, sending}

var I_lyr_xslvar { li in XSLVAR, p in PHASES} >= -2.0 <= 2.0;
		# Imag current of a xslvar's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lys_xslvar { li in XSLVAR, p in PHASES} >= -2.0 <= 2.0; 
# Real current of a xslvar's
        # e end connected to bus n phase p. e in {receiving, sending}

# var I_lr2 { li in XSLVAR, p in PHASES} =  I_lxr [li, p]^2 + I_lyr [li, p]^2 ;
var I_lrmag_xslvar { li in XSLVAR, p in PHASES} 
	=  sqrt (I_lxr_xslvar [li, p]^2 + I_lyr_xslvar [li, p]^2 );
var I_lrangle_xslvar { li in XSLVAR, p in PHASES} 
	=  atan2(I_lyr_xslvar [li, p], I_lxr_xslvar [li, p])*180/pi ;

var xslvar_pr { li in XSLVAR, p in PHASES} =
     Ex[xslvar_rbus[li],p]*I_lxr_xslvar[li, p] 
	+ Ey[xslvar_rbus[li],p]*I_lyr_xslvar[li, p] ; 

var xslvar_ps { li in XSLVAR, p in PHASES} =
     Ex[xslvar_sbus[li],p]*I_lxs_xslvar[li, p] 
	+ Ey[xslvar_sbus[li],p]*I_lys_xslvar[li, p] ; 

var xslvar_loss { li in XSLVAR, p in PHASES} = 
	xslvar_ps[li, p] - xslvar_pr[li, p];
	
var all_xslvar_loss = sum { li in XSLVAR, p in PHASES}
                     xslvar_loss[li, p];

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##  XSLVAR' real voltage constraints	
subject to xslvar_real_volt{ li in XSLVAR, p in PHASES}:
    sum {pa in PHASES} ( xslvar_Ax[li, p, pa] *Ex[xslvar_rbus[li],pa]
                                    - xslvar_Ay[li, p, pa]*Ey[xslvar_rbus[li], pa]
                                    + xslvar_Bx[li, p, pa]*I_lxr_xslvar[li, pa]
                                    - xslvar_By[li, p, pa]*I_lyr_xslvar[li, pa] )
    =Ex[xslvar_sbus[li],p];
## XSLVAR' imag voltage constraints
subject to xslvar_imag_volt{ li in XSLVAR, p in PHASES}:
    sum {pa in PHASES} ( xslvar_Ax[li, p, pa] *Ey[xslvar_rbus[li],pa] 
	+  xslvar_Ay[li, p, pa]*Ex[xslvar_rbus[li], pa] 
    +  xslvar_Bx[li, p, pa]*I_lyr_xslvar[li, pa] 
    +  xslvar_By[li, p, pa]*I_lxr_xslvar[li, pa] )
        =Ey[xslvar_sbus[li],p];	

subject to xslvar_uniform12 {xsl in XSLVAR}:
	xslvar_By[xsl, 1, 1] = xslvar_By[xsl, 2, 2];
subject to xslvar_unifor13 {xsl in XSLVAR}:
	xslvar_By[xsl, 1, 1] = xslvar_By[xsl, 3, 3];

subject to xslvar_By_bin_constraint {xsl in XSLVAR, p in PHASES}:
	xslvar_By[xsl, p, p] = xslvar_By_lower_bound[xsl,p] + xslvar_By_bin [xsl,p]*0.01;

## xslvar' imag current constraints
subject to xslvar_imag_current{ li in XSLVAR, p in PHASES}:
    sum {pa in PHASES} ( xslvar_C[li, p, pa] *Ey[xslvar_rbus[li],pa]
                                    + xslvar_D[li, p, pa]*I_lyr_xslvar[li, pa] )
        =I_lys_xslvar[li,p];
## Lines' real current constraints
subject to xslvar_real_current{ li in XSLVAR, p in PHASES}:
    sum {pa in PHASES} ( xslvar_C[li, p, pa] *Ex[xslvar_rbus[li],pa]
                                    + xslvar_D[li, p, pa]*I_lxr_xslvar[li, pa] )
        =I_lxs_xslvar[li,p];			

## Real current balance equalities
subject to current_balance_real_xslvar {n in BUS,p in PHASES,
       xsl in XSLVAR: xslvar_rbus[xsl] = n}:
    I_lxr_xslvar[xsl,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lxs[lin,p]
	# from loads connected to this node

	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Lx[L,p] + I_LxD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Lx_Imp[L,p] + I_LxD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Lx_Cur[L,p] + I_LxD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Lx_dg[L,p] + I_LxD_dg[L,p])

	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Ix[s,p];		
	# from other XSLVAR
	;

		
## imag current balance equalities
subject to current_balance_imag_xslvar {n in BUS, p in PHASES, xsl in XSLVAR: 
     xslvar_rbus[xsl] = n}:
    I_lyr_xslvar[xsl,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lys[lin,p] 
	# from loads connected to this node
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Ly[L,p] + I_LyD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Ly_Imp[L,p] + I_LyD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Ly_Cur[L,p] + I_LyD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Ly_dg[L,p] + I_LyD_dg[L,p])

	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Iy[s,p];	

	

################################
# End of XSLVAR Model
################################
################################
# OLTC Model requires gen model
################################
# include OLTC.mod;
# OLTC Model
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
set OLTC_TAPS := {-0.1    , -0.09375, -0.0875 , -0.08125, -0.075  , -0.06875,
       -0.0625 , -0.05625, -0.05   , -0.04375, -0.0375 , -0.03125,
       -0.025  , -0.01875, -0.0125 , -0.00625,  0.     ,  0.00625,
        0.0125 ,  0.01875,  0.025  ,  0.03125,  0.0375 ,  0.04375,
        0.05   ,  0.05625,  0.0625 ,  0.06875,  0.075  ,  0.08125,
        0.0875 ,  0.09375,  0.1    };
set OLTC_TAPS_BINARY :={0,1};
param oltc_lower_bound {o in OLTC, p in PHASES};
param oltc_name{OLTC} symbolic;
param oltc_type{OLTC}; # types:
                # 1. default
param oltc_number{OLTC}; # number of a oltc
param oltc_Ax{OLTC, PHASES, PHASES}; # real part of A matrix

param oltc_Ay{OLTC, PHASES, PHASES}; # imaginary part of A matrix
param oltc_Bx{OLTC, PHASES, PHASES}; # real part of B matrix
param oltc_By{OLTC, PHASES, PHASES}; # imaginary part of B matrix
param oltc_C{OLTC, PHASES, PHASES}; # real part of C matrix
param oltc_Bmag {li in OLTC, p1 in PHASES, p2 in PHASES} := sqrt ( oltc_Bx[li, p1, p2]^2 + oltc_By[li, p1, p2]^2) ;
param oltc_R {li in OLTC, p1 in PHASES, p2 in PHASES: oltc_Bmag[li, p1, p2] != 0} := oltc_Bx[li, p1, p2]/oltc_Bmag[li, p1, p2] ;
##param oltc_Cx{LINES, PHASES, PHASES}; # real part of C matrix
##param oltc_Cy{LINES, PHASES, PHASES}; # imag part of C matrix
param oltc_D{OLTC, PHASES, PHASES}; # real part of D matrix
##param oltc_Dx{LINES, PHASES, PHASES}; # real part of D matrix
##param oltc_Dy{LINES, PHASES, PHASES}; # imag part of D matrix
param oltc_sbus{li in OLTC}; # a oltc's sending bus
param oltc_rbus{li in OLTC}; # a oltc's receiving bus


# variable for OLTC
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
var I_lxr_oltc { li in OLTC, p in PHASES} >= -2.0 <= 2.0; # Real current of a oltc's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lxs_oltc { li in OLTC, p in PHASES} >= -2.0 <= 2.0; # Real current of a oltc's
        # e end connected to bus n phase p. e in {receiving, sending}

var I_lyr_oltc { li in OLTC, p in PHASES} >= -2.0 <= 2.0;
		# Imag current of a oltc's
        # e end connected to bus n phase p. e in {receiving, sending}
var I_lys_oltc { li in OLTC, p in PHASES} >= -2.0 <= 2.0; # Real current of a oltc's
        # e end connected to bus n phase p. e in {receiving, sending}

var I_lr2_oltc { li in OLTC, p in PHASES} =  I_lxr_oltc [li, p]^2 + I_lyr_oltc [li, p]^2 ;
var I_lrmag_oltc { li in OLTC, p in PHASES} =  sqrt (I_lxr_oltc [li, p]^2 + I_lyr_oltc [li, p]^2 );
var I_lrangle_oltc { li in OLTC, p in PHASES} =  atan2(I_lyr_oltc [li, p], I_lxr_oltc [li, p])*180/pi ;

# var oltc_loss { li in LINES, p in PHASES} = 
	# sum {pa in PHASES : oltc_Bmag[li, p, pa] != 0} 
	# ( oltc_Bx[li, p, pa] * I_lr2[li, pa]/oltc_Bmag[li, p, pa] ) ;


## losses of oltc 
var oltc_loss_oltc { li in OLTC, p in PHASES} = 
	sum {pa in PHASES : oltc_Bmag[li, p, pa] != 0 and p=pa} 
	( oltc_Bx[li, p, pa] * I_lr2[li, pa]) ;


var oltc_pr { li in OLTC, p in PHASES} =
     Ex[oltc_rbus[li],p]*I_lxr_oltc[li, p] 
	+ Ey[oltc_rbus[li],p]*I_lyr_oltc[li, p] ; 

var oltc_taps_bin {t in OLTC, p in PHASES}  in OLTC_TAPS_BINARY;#>= -0.1 <= 0.1;
var oltc_taps {t in OLTC, p in PHASES} =  oltc_lower_bound[t, p] + oltc_taps_bin[t, p]*0.00625;

# ##  OLTCs' real voltage constraints	
subject to oltc_real_volt{ li in OLTC, p in PHASES}:
    sum {pa in PHASES: p=pa} ( Ex[oltc_rbus[li],pa]/(oltc_Ax[li, p, pa] + oltc_taps[li, pa]) 
	                                       
            - oltc_Ay[li, p, pa]*Ey[oltc_rbus[li], pa]
            + oltc_Bx[li, p, pa]*I_lxr_oltc[li, pa]
            - oltc_By[li, p, pa]*I_lyr_oltc[li, pa] )
    =Ex[oltc_sbus[li],p];
## 6. OLTCs' imag voltage constraints. oltc_A is a real only matrix. Hence I do 
## not add oltc_taps part as in oltc_Ax part.
subject to oltc_imag_volt { li in OLTC, p in PHASES}:
    sum {pa in PHASES: p=pa} ( Ey[oltc_rbus[li],pa]/(oltc_Ax[li, p, pa] + oltc_taps[li, pa])
	+  oltc_Ay[li, p, pa]*Ex[oltc_rbus[li], pa] 
    +  oltc_Bx[li, p, pa]*I_lyr_oltc[li, pa] 
    +  oltc_By[li, p, pa]*I_lxr_oltc[li, pa] )
        =Ey[oltc_sbus[li],p];	
		
## OLTCs' imag current constraints
subject to oltc_imag_current{ li in OLTC, p in PHASES}:
    sum {pa in PHASES} ( oltc_C[li, p, pa] *Ey[oltc_rbus[li],pa]
                                    + (oltc_D[li, p, pa]+oltc_taps[li,pa])*I_lyr_oltc[li, pa] )
        =I_lys_oltc[li,p];
## OLTCs' real current constraints
subject to oltc_real_current{ li in OLTC, p in PHASES}:
    sum {pa in PHASES: p= pa} ( oltc_C[li, p, pa] *Ex[oltc_rbus[li],pa]
                                    + (oltc_D[li, p, pa]+ oltc_taps[li,pa])*I_lxr_oltc[li, pa] )
        =I_lxs_oltc[li,p];



		
## Real current balance equalities
subject to current_balance_real_oltc {n in BUS, p in PHASES, ol in OLTC: 
     oltc_rbus[ol] = n}:
    I_lxr_oltc[ol,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lxs[lin,p]
	# from loads connected to this node
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Lx[L,p] + I_LxD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Lx_Imp[L,p] + I_LxD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Lx_Cur[L,p] + I_LxD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Lx_dg[L,p] + I_LxD_dg[L,p])
	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Ix[s,p]
	;
	
## imag current balance equalities
subject to current_balance_imag_oltc {n in BUS, p in PHASES, ol in OLTC: 
     oltc_rbus[ol] = n}:
    I_lyr_oltc[ol,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lys[lin,p] 
	# from loads connected to this node
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Ly[L,p] + I_LyD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Ly_Imp[L,p] + I_LyD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Ly_Cur[L,p] + I_LyD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Ly_dg[L,p] + I_LyD_dg[L,p])
	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Iy[s,p];	

# the following to constraints ensure uniform operation of OLTC.
# if non uniform operation is needed, the they can be deactivated
# subject to OLTC_uniform1 {ol in OLTC}:
	# oltc_taps[ol, 1] - oltc_taps[ol, 2] = 0;
# subject to OLTC_uniform2 {ol in OLTC}:
	# oltc_taps[ol, 1] - oltc_taps[ol, 3] = 0;
################################
# End of OLTC Model
################################



######################
# objective function #
######################

 var voltage_error_squared = 
		# sum {n in BUS, p in PHASES} 
      sum {n in BUS, p in PHASES:  node_type[n, p] = 2} 
	     ( 1.- sqrt(Ex[n,p]^2 + Ey[n,p]^2 +1e-12) )^2;

minimize voltage_deviation:
    voltage_error_squared;
# minimize voltage_error_squared:
    # sum {n in BUS, p in PHASES:  node_type[n, p] = 2} 
	# # sum {n in BUS, p in PHASES} 
	     # ( 1.- sqrt(Ex[n,p]^2 + Ey[n,p]^2));#+1e-012));
minimize line_loss_obj:
	# sum {li in LINES, p in PHASES}	 line_loss[li,p];
	# system_losses;
	all_line_loss;

minimize Psub_obj:
	Psub;

minimize Qsub_obj:
	Qsub;
	
minimize PandQsub_obj:
	Psub+Qsub;

minimize all_obj :
   all_line_loss + Psub + Qsub  + voltage_error_squared;
################################
# End of Objective Functions
################################


###############
# Constraints #
###############
# Equality constraints
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Lines' imag current constraints
subject to line_imag_current{ li in LINES, p in PHASES}:
    sum {pa in PHASES} ( line_C[li, p, pa] *Ey[line_rbus[li],pa]
                                    + line_D[li, p, pa]*I_lyr[li, pa] )
        =I_lys[li,p];
## Lines' real current constraints
subject to line_real_current{ li in LINES, p in PHASES}:
    sum {pa in PHASES} ( line_C[li, p, pa] *Ex[line_rbus[li],pa]
                                    + line_D[li, p, pa]*I_lxr[li, pa] )
        =I_lxs[li,p];

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~		
## Real current balance equalities
subject to current_balance_real {n in BUS, p in PHASES, li in LINES : 
     line_rbus[li] = n}:
    I_lxr[li,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lxs[lin,p]
	# from loads connected to this node
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Lx[L,p] + I_LxD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Lx_Imp[L,p] + I_LxD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Lx_Cur[L,p] + I_LxD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Lx_dg[L,p] + I_LxD_dg[L,p])
	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Ix[s,p]
	# from other OLTC connected to this node. 
	+ sum {ol in OLTC: oltc_sbus[ol] = n} I_lxs_oltc[ol,p]
	# from XSLVAR connected to this OLTC.
	+ sum {xsl in XSLVAR: xslvar_sbus[xsl] = n} I_lxs_xslvar[xsl,p]
	;
	
	
## imag current balance equalities

subject to current_balance_imag {n in BUS, p in PHASES, li in LINES : 
     line_rbus[li] = n}:
    I_lyr[li,p] =
	# from lines below this node
	  sum {lin in LINES: line_sbus[lin] = n}I_lys[lin,p] 
	# from loads connected to this node
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=1}  (I_Ly[L,p] + I_LyD[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=2}  (I_Ly_Imp[L,p] + I_LyD_Imp[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=3}  (I_Ly_Cur[L,p] + I_LyD_Cur[L,p])
	+ sum {L in LOAD : load_bus[L] = n and load_type[L]=4}  (I_Ly_dg[L,p] + I_LyD_dg[L,p])
	# from generator connected to this node
	+ sum {s in GEN: gen_bus[s] = n} gen_Iy[s,p]
    # from other OLTC connected this node.
	+ sum {ol in OLTC: oltc_sbus[ol] = n} I_lys_oltc[ol,p]
	# from XSLVAR connected to this OLTC.
	+ sum {xsl in XSLVAR: xslvar_sbus[xsl] = n} I_lys_xslvar[xsl,p]	
	;	
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ##  Lines' real voltage constraints	
subject to line_real_volt{ li in LINES, p in PHASES}:
    sum {pa in PHASES} ( line_Ax[li, p, pa] *Ex[line_rbus[li],pa]
                                    - line_Ay[li, p, pa]*Ey[line_rbus[li], pa]
                                    + line_Bx[li, p, pa]*I_lxr[li, pa]
                                    - line_By[li, p, pa]*I_lyr[li, pa] )
    =Ex[line_sbus[li],p];
## 6. Lines' imag voltage constraints
subject to line_imag_volt{ li in LINES, p in PHASES}:
    sum {pa in PHASES} ( line_Ax[li, p, pa] *Ey[line_rbus[li],pa] 
	+  line_Ay[li, p, pa]*Ex[line_rbus[li], pa] 
    +  line_Bx[li, p, pa]*I_lyr[li, pa] 
    +  line_By[li, p, pa]*I_lxr[li, pa] )
        =Ey[line_sbus[li],p];	

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~		
## SVC's XSL constraints 
# ##  Lines' real voltage constraints	for XSL ===> static version
subject to line_real_volt_xsl{ li in LINES, p in PHASES: line_type[li] = 4  }:
    sum {pa in PHASES : p=pa} ( line_Ax[li, p, pa] *Ex[line_rbus[li],pa]
        - line_Ay[li, p, pa]*Ey[line_rbus[li], pa]
        + line_Bx[li, p, pa]*I_lxr[li, pa]
        - line_By[li, p, pa]*I_lyr[li, pa]
	#	+ line_xsl[li]*I_lyr[li, pa] 
					   )
    =Ex[line_sbus[li],p];
	
## 6. Lines' imag voltage constraints XSL ===> static version
subject to line_imag_volt_xsl{ li in LINES, p in PHASES : line_type[li] = 4 }:
    sum {pa in PHASES: p=pa} ( line_Ax[li, p, pa] *Ey[line_rbus[li],pa] 
	+  line_Ay[li, p, pa]*Ex[line_rbus[li], pa] 
    +  line_Bx[li, p, pa]*I_lyr[li, pa] 
    +  line_By[li, p, pa]*I_lxr[li, pa]
	#+  line_xsl[li]*I_lxr[li, pa]  
	)
        =Ey[line_sbus[li],p];	

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# GEN's constraints
subject to P_balance {s in GEN, p in PHASES}:
	  Ex[gen_bus[s],p]*gen_Ix[s, p] 
	+ Ey[gen_bus[s],p]*gen_Iy[s, p] = 0;

# subject to PV_V {s in GEN, p in PHASES}:
	 # Emag[gen_bus[s],p] = gen_vref[s];
	 
# positive sequence voltag magnitude equality for a PV bus.
subject to PV_V {s in GEN}:
	 E1mag[s] = gen_vref[s];

# These can be used if we requires each phase
# generates the same amount of reactive power.
# however, it can not be used at the same time
# with the voltage equality constraints because
# the unbalance of lines and loads.
subject to gen_q12 {s in GEN}:
	 gen_q[s,1] = gen_q[s,2] ;
subject to gen_q13 {s in GEN}:
	 gen_q[s,1] = gen_q[s,3] ;	 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# # Inequality constraints	
# subject to Ex_c1 {n in BUS: node_type[n, 1] = 2}:
	# 0.0 <= Ex[n, 1] <= 1.0;

# subject to Ex_c2 {n in BUS: node_type[n, 2] = 2}:
	# -1.0 <= Ex[n, 2] <= 0.0;

# subject to Ex_c3 {n in BUS: node_type[n, 3] = 2}:
	# -1.0 <= Ex[n, 3] <= 0.0;	

# subject to Ey_c1 {n in BUS: node_type[n, 1] = 2}:
	# -1.0 <= Ey[n, 1 ] <= 1.0;
	
# subject to Ey_c2 {n in BUS : node_type[n, 2] = 2}:
	# -1.0 <= Ey[n, 2 ] <= 0.0;
	
# subject to Ey_c3 {n in BUS: node_type[n, 3] = 2}:
	# 0 <= Ey[n, 3 ] <= 1.0;
	
subject to Emag_c {n in BUS, p in PHASES: node_type[n, p] = 2}:
# subject to Emag_c {n in BUS, p in PHASES}:
	# 0.9025 <= Ex[n, p]^2+Ey[n, p]^2 <= 1.1025;
	0.95 <= sqrt(Ex[n, p]^2+Ey[n, p]^2) <= 1.05;

subject to gen_qc {g in GEN, p in PHASES}:
	-0.12 <= gen_q[g, p] <= 0.05;
# subject to line_xslc {li in LINES : line_type[li]=4}:
	# 0.0 <= line_xsl[li] <= 0.05;
# generator voltage reference limits. This must be an appropriate value
# considering the fact that the voltage at source bus is fixed. I think
# it should be it such a way that it does not try to increase the fixed
# source bus.
subject to gen_vref_c {s in GEN}:
    0.95<= gen_vref[s] <= 1.05;

subject to xslvar_By_c {x in XSLVAR, p in PHASES}:
	 0.0<= xslvar_By[x, p, p] <= 0.05;

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
################################
# End of COSTRAINTS
################################


display '===========================================================';
objective line_loss_obj;