From 2f8c5e47435fb36305f31e5da7daec42477b432f Mon Sep 17 00:00:00 2001
From: ckhanpour3 <ckhanpour3@gatech.edu>
Date: Mon, 4 Nov 2024 16:34:37 -0500
Subject: [PATCH] AC Power Flow Implementation

First attempt on AC Power Flow constraints in gtep_model.py

Added max/min reactive power constraints for thermal generators.

Added "ACP" config option for supported flows in config_options.py
---
 gtep/config_options.py |  2 ++
 gtep/gtep_model.py     | 76 +++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 77 insertions(+), 1 deletion(-)

diff --git a/gtep/config_options.py b/gtep/config_options.py
index 5bdea95..2158b32 100644
--- a/gtep/config_options.py
+++ b/gtep/config_options.py
@@ -14,6 +14,8 @@
 _supported_flows = {
     "DC": ("gtep.dcopf", "DC power flow approximation"),
     "CP": ("gtep.cp", "Copper plate power flow approximation"),
+    "ACP": ("gtep.acp", "AC power flow in polar formulation"),
+    
 }
 
 
diff --git a/gtep/gtep_model.py b/gtep/gtep_model.py
index ce87db0..4f44b40 100644
--- a/gtep/gtep_model.py
+++ b/gtep/gtep_model.py
@@ -636,6 +636,60 @@ def dc_power_flow(disj):
             return b.powerFlow[branch] == (-1 / reactance) * (
                 disj.busAngle[tb] - disj.busAngle[fb] + shift
             )
+        
+        if m.config["flow_model"] == "ACP":
+            fb = m.transmission[branch]["from_bus"]
+            tb = m.transmission[branch]["to_bus"]
+            resistance = m.md.data["elements"]["branch"][branch].get("resistance", 0.0)
+            reactance = m.md.data["elements"]["branch"][branch].get("reactance", 1e-6) 
+
+            # Transformer tap ratio and phase shift 
+            if m.md.data["elements"]["branch"][branch]["branch_type"] == "transformer":
+                reactance *= m.md.data["elements"]["branch"][branch][
+                    "transformer_tap_ratio"
+                ]
+                phase_shift = m.md.data["elements"]["branch"][branch][
+                    "transformer_phase_shift"
+                ]
+            else:
+                phase_shift = 0
+
+            admittance = 1 / complex(resistance, reactance)
+            G = admittance.real
+            B = admittance.imag
+            
+            # Define voltage magnitude variables for from and to buses
+            @disj.Var()
+            def voltage_from(disj):
+                return disj.parent_block().voltage[fb]
+
+            @disj.Var()
+            def voltage_to(disj):
+                return disj.parent_block().voltage[tb]
+
+
+            # Active Power Flow Constraint
+            @disj.Constraint()
+            def ac_power_flow_p(disj):
+                return m.P_flow[branch] == (
+                    m.voltage[fb]**2 * G
+                    - m.voltage[fb] * m.voltage[tb] * (
+                        G * math.cos(m.theta[fb] - m.theta[tb] + phase_shift)
+                        + B * math.sin(m.theta[fb] - m.theta[tb] + phase_shift)
+                    )
+                )
+            
+            # Reactive Power Flow Constraint
+            @disj.Constraint()
+            def ac_power_flow_q(disj):
+                return m.Q_flow[branch] == (
+                    -m.voltage[fb]**2 * B
+                    - m.voltage[fb] * m.voltage[tb] * (
+                        G * math.sin(m.theta[fb] - m.theta[tb] + phase_shift)
+                        - B * math.cos(m.theta[fb] - m.theta[tb] + phase_shift)
+                    )
+                )
+
 
     @b.Disjunct(m.transmission)
     def branchNotInUse(disj, branch):
@@ -755,7 +809,7 @@ def capacity_factor(b, renewableGen):
             b.renewableGeneration[renewableGen] + b.renewableCurtailment[renewableGen]
             == m.renewableCapacity[renewableGen]
         )
-
+    
 
     ## TODO: (@jkskolf) add renewableExtended to this and anywhere else
     @b.Constraint(m.renewableGenerators)
@@ -1568,6 +1622,26 @@ def model_data_references(m):
         for thermalGen in m.thermalGenerators
     }
 
+    # Maximum reactive power output of each thermal generator
+    m.thermalReactiveCapacity = {
+        thermalGen: m.md.data["elements"]["generator"][thermalGen].get("q_max", 0)
+        for thermalGen in m.thermalGenerators
+    }
+
+    # Minimum reactive power output of each thermal generator
+    m.thermalReactiveMin = {
+        thermalGen: m.md.data["elements"]["generator"][thermalGen].get("q_max", 0)
+        for thermalGen in m.thermalGenerators
+    }  
+
+    # Handles case where reactive power limits of generator is not provided
+    for gen in m.thermalGenerators:
+        if gen not in m.thermalReactiveCapacity:
+            m.thermalReactiveCapacity[gen] = 0
+        if gen not in m.thermalReactiveMin:
+            m.thermalReactiveMin[gen] = 0
+
+
     # Maximum output of each renewable generator
     m.renewableCapacity = {
         renewableGen: (