Combined the extend.jl and power_flow.jl files

src/REopt.jl

@@ -197,7 +197,6 @@ include("outagesim/outage_simulator.jl")
 include("outagesim/backup_reliability.jl")
 
 include("lindistflow/power_flow.jl")
-include("lindistflow/extend.jl")
 
 include("mpc/results.jl")
 include("mpc/model.jl")

src/constraints/electric_utility_constraints.jl

@@ -3,15 +3,15 @@ function add_export_constraints(m, p; _n="")
 
     ##Constraint (8e): Production export and curtailment no greater than production
     if string(p.s.site.node) != p.s.settings.facilitymeter_node
-        print("\n Updated***Test 4,  Adding constraint 8e to node $(p.s.site.node)")
+        print("\n Updated***4000 timesteps,  Adding constraint 8e to node $(p.s.site.node)")
         @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
             p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] 
             >= sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t]) +
             m[Symbol("dvCurtail"*_n)][t, ts]
         )
     else
         print("\n Not adding constraint 8e to the facility meter node, node $(p.s.site.node)")
-        TempVector = collect(201:8760) # implemented only for initial testing; TODO: remove this after testing
+        TempVector = collect(4001:8760) # implemented only for initial testing; TODO: remove this after testing
         @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid], m[Symbol("dvRatedProduction"*_n)][t,ts] == 0)
         @constraint(m, [t in p.techs.elec, ts in TempVector], sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t]) == 0)
         @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],  m[Symbol("dvCurtail"*_n)][t, ts] == 0)

src/lindistflow/power_flow.jl

@@ -1,4 +1,6 @@
-
+#=
+Acknowledgement: This code is based on code from the LinDistFlow.jl package
+=#
 
 """
 mutable struct PowerFlowInputs
@@ -183,6 +185,152 @@ PowerFlowInputs(
 )
 end
 
+function build_power_flow!(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1};
+    make_import_export_complementary::Bool=true)
+
+    power_flow_add_variables(m, p)
+    add_expressions(m, ps)
+    constrain_power_balance(m, p)
+    constrain_substation_voltage(m, p)
+    create_line_variables(m, p)
+    constrain_loads(m, p, ps)
+    # Note: the constrain_KVL(m, p) function is called in the microgrid.jl file
+    if make_import_export_complementary
+        add_complementary_constraints(m, ps)
+    end
+end
+
+#2. Add power flow constraints to m, 
+#    - set Pⱼ's = -1 * (dvGridPurchase_j - dvProductionToGrid_j)
+function add_expressions(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1})
+for p in ps
+    _n = string("_", p.s.site.node)
+    if string(p.s.site.node) != p.s.settings.facilitymeter_node
+        print("\n Applying total export equation for node $(p.s.site.node)")
+        m[Symbol("TotalExport"*_n)] = @expression(m, [ts in p.time_steps],
+            sum(
+                m[Symbol("dvProductionToGrid"*_n)][t,u,ts] 
+                for t in p.techs.elec, u in p.export_bins_by_tech[t]
+            )
+            + sum(m[Symbol("dvStorageToGrid"*_n)][b,ts] for b in p.s.storage.types.all )# added this line to include battery export in the total export
+        )
+    else
+        # set the total node export to 0 for the facility meter grid, because that node has no techs
+            # also, for that node, the dvProductionToGrid is used for the grid export benefits and set to the powerflow of the substation line when flow on that line is negative
+        print("\n Setting the total export for node $(p.s.site.node) to zero. This node is the facility meter node.")
+
+        dv = "TotalExport_"*p.s.settings.facilitymeter_node
+        m[Symbol(dv)] = @variable(m, [p.time_steps], base_name=dv, lower_bound =0)
+
+        @constraint(m, [t in p.time_steps], m[Symbol("TotalExport_"*p.s.settings.facilitymeter_node)][t] == 0)
+
+    end
+end
+end
+
+# TODO add complementary constraint to UL for dvProductionToGrid_ and dvGridPurchase_ (don't want it in LL s.t. it stays linear)
+function add_complementary_constraints(m::JuMP.AbstractModel, ps::Array{REoptInputs{Scenario}, 1})
+for p in ps
+    _n = string("_", p.s.site.node)
+        print("\n Adding the complementary constraint to node $(p.s.site.node)")
+        for (i, e) in zip(m[Symbol("dvGridPurchase"*_n)], m[Symbol("TotalExport"*_n)])
+            @constraint(m,
+                [i, e] in MOI.SOS1([1.0, 2.0])
+            )
+        end
+end
+end
+
+
+function constrain_loads(m::JuMP.AbstractModel, p::PowerFlowInputs, ps::Array{REoptInputs{Scenario}, 1})
+
+Pⱼ = m[:Pⱼ]
+Qⱼ = m[:Qⱼ]
+
+reopt_nodes =  [p.s.site.node for p in ps] 
+
+# Note: positive values are injections
+
+for j in p.busses
+    if j in keys(p.Pload)
+        #print("\n Debugging: the j variable is: $(j) ")
+        if parse(Int, j) in reopt_nodes
+            if j != "15" 
+                #print("\n Debuggin: j is not 15")               
+                @constraint(m, [t in 1:p.Ntimesteps],
+                    Pⱼ[j,t] == 1e3/p.Sbase * (  # 1e3 b/c REopt values in kW
+                        m[Symbol("TotalExport_" * j)][t]
+                        - m[Symbol("dvGridPurchase_" * j)][t]
+                    )
+                )
+            else
+                print("\n j is 15 and the j variable is: $(j)")
+                @constraint(m, [t in 1:p.Ntimesteps], Pⱼ["15",t] == 0)
+            end
+        else
+            # This constraint is for the power balance for nodes that have loads defined through ldf, but not in the REopt inputs
+            @constraint(m, [t in 1:p.Ntimesteps],
+                Pⱼ[j,t] == -p.Pload[j][t]
+            )
+        end
+    elseif j != p.substation_bus
+        @constraint(m, [t in 1:p.Ntimesteps],
+            Pⱼ[j,t] == 0
+        )
+    end
+
+    if j in keys(p.Qload)
+        if parse(Int, j) in reopt_nodes
+            if j != "15"
+                #print("\n Debuggin: (reactive power) j is not 15") 
+                @constraint(m, [t in 1:p.Ntimesteps],
+                    Qⱼ[j,t] == 1e3/p.Sbase * p.pf * (  # 1e3 b/c REopt values in kW
+                    m[Symbol("TotalExport_" * j)][t]
+                        - m[Symbol("dvGridPurchase_" * j)][t]
+                    )
+                )
+            else
+                print("\n (reactive power) j is 15 and the j variable is: $(j)")
+                @constraint(m, [t in 1:p.Ntimesteps], Qⱼ["15",t] == 0)
+            end
+        else
+            @constraint(m, [t in 1:p.Ntimesteps],
+                Qⱼ[j,t] == -p.Qload[j][t]
+            )
+        end
+    elseif j != p.substation_bus
+        @constraint(m, [t in 1:p.Ntimesteps],
+            Qⱼ[j,t] == 0
+        )
+    end
+end
+p.Nequality_cons += 2 * (p.Nnodes - 1) * p.Ntimesteps
+
+# Constrain loads on the transformers
+P = m[:Pᵢⱼ] 
+Q = m[:Qᵢⱼ]
+all_transformers = []
+# Define decision variable for the transformer maximum kVa
+for i in keys(p.transformers)
+    if p.transformers[i]["Transformer Side"] == "downstream"
+        push!(all_transformers, i)
+    end
+end
+print("\n The all_transformers variable is: ")
+print(all_transformers)
+@variable(m, transformer_max_kva[all_transformers] >= 0)
+
+# Apply the transformer max kva to the constraints
+for i in keys(p.transformers)
+    if p.transformers[i]["Transformer Side"] == "downstream"           
+            DirectlyUpstreamNode = i_to_j(i, p)  
+            transformer_internal_line = string(DirectlyUpstreamNode[1])*"-"*string(i)
+            @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] <= ((m[:transformer_max_kva][i]*1000)/p.Sbase))
+            @constraint(m, [T in 1:p.Ntimesteps], P[transformer_internal_line, T] + Q[transformer_internal_line, T] >= -((m[:transformer_max_kva][i]*1000)/p.Sbase))
+    end
+end
+end
+
 
 function power_flow_add_variables(m, p::PowerFlowInputs)
     T = 1:p.Ntimesteps