Draft Pull Request: Multinode powerflow updates

Project.toml

@@ -11,35 +11,35 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
-LinDistFlow = "bf674bac-ffe4-48d3-9f32-72124ffa9ede"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b"
 
 [compat]
 ArchGDAL = "0.8, 0.9"
-CoolProp = "0.1"
 CSV = "0.10"
+CoolProp = "0.1"
 DataFrames = "1.4"
+Dates = "1"
+DelimitedFiles = "1"
 HTTP = "0.8, 0.9, 1"
 JLD = "0.13"
 JSON = "0.21"
 JuMP = "0.21, 0.22, 0.23, 1"
-LinDistFlow = "0.1, 0.2"
+LinearAlgebra = "1"
+Logging = "1"
 MathOptInterface = "0.9, 0.10, 1"
 Requires = "1.3"
 Roots = "1.3, 2"
-TestEnv = "1.7, 1.8, 1.9"
-Dates = "1"
-DelimitedFiles = "1"
-LinearAlgebra = "1"
-Logging = "1"
 Statistics = "1"
+TestEnv = "1.7, 1.8, 1.9"
 julia = "1.4, 1.5, 1.6, 1.7, 1.8"

src/REopt.jl

@@ -11,7 +11,6 @@ export
     simulate_outages,
     add_variables!,
     add_objective!,
-    LinDistFlow,
     MPCScenario,
     MPCInputs,
     run_mpc,
@@ -30,9 +29,6 @@ export
 
 import HTTP
 import JSON
-using LinDistFlow  # required to export LinDistFlow
-import LinDistFlow 
-const LDF = LinDistFlow
 using JuMP
 using JuMP.Containers: DenseAxisArray
 using Logging
@@ -52,9 +48,14 @@ using CoolProp
 using LinearAlgebra
 using CSV
 using DataFrames
+using Plots
+import Plots
+using HiGHS
 
 function __init__()
     @require GhpGhx="7ce85f02-24a8-4d69-a3f0-14b5daa7d30c" println("using GhpGhx module in REopt")
+    @require PowerModelsDistribution="d7431456-977f-11e9-2de3-97ff7677985e" println("using the PowerModelsDistribution package") # This will load whichever feature branch you have installed in the coding environment
+    @require Xpress="9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054" println("using the Xpress package")
 end
 
 const EXISTING_BOILER_EFFICIENCY = 0.8
@@ -195,11 +196,15 @@ include("results/electric_heater.jl")
 include("results/ashp.jl")
 
 include("core/reopt.jl")
-include("core/reopt_multinode.jl")
 include("outagesim/outage_simulator.jl")
 include("outagesim/backup_reliability.jl")
 
-include("lindistflow/extend.jl")
+include("multinode/reopt_multinode.jl")
+include("multinode/multinode.jl")
+include("multinode/multinode_outage_simulator.jl")
+include("multinode/multinode_results_processing.jl")
+include("multinode/multinode_inputs.jl")
+include("multinode/multinode_utils.jl")
 
 include("mpc/results.jl")
 include("mpc/model.jl")

src/constraints/electric_utility_constraints.jl

@@ -1,12 +1,17 @@
 # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE.
 function add_export_constraints(m, p; _n="")
-
     ##Constraint (8e): Production export and curtailment no greater than production
-    @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]
-    )
+    if (_n == "") || (p.s.settings.facilitymeter_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
+        # Don't add constraint 8e to the facility meter node
+        @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 p.time_steps_with_grid],  m[Symbol("dvCurtail"*_n)][t, ts] == 0)
+    end
 
     binNEM = 0
     binWHL = 0
@@ -19,12 +24,16 @@ function add_export_constraints(m, p; _n="")
     if !isempty(NEM_techs)
         # Constraint (9c): Net metering only -- can't sell more than you purchase
         # hours_per_time_step is cancelled on both sides, but used for unit consistency (convert power to energy)
-        @constraint(m,
-            p.hours_per_time_step * sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
-            for t in NEM_techs, ts in p.time_steps)
-            <= p.hours_per_time_step * sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
-                for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)
-        )
+        if (_n == "") || (p.s.settings.facilitymeter_node != "")
+            @constraint(m,
+                p.hours_per_time_step * sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
+                for t in NEM_techs, ts in p.time_steps)
+                <= p.hours_per_time_step * sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
+                    for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)
+            )
+        else
+            # Don't add constraint 9c to the facility meter node
+        end
 
         if p.s.electric_utility.net_metering_limit_kw == p.s.electric_utility.interconnection_limit_kw && isempty(WHL_techs)
             # no need for binNEM nor binWHL
@@ -134,9 +143,17 @@ function add_export_constraints(m, p; _n="")
 
         if typeof(binNEM) <: Real  # no need for wholesale binary
             binWHL = 1
+            for b in p.s.storage.types.all
+                println(b)
+            end
             WHL_benefit = @expression(m, p.pwf_e * p.hours_per_time_step *
-                sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
-                        for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                  (
+                    sum(sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvStorageToGrid"*_n)][b, ts] 
+                        for b in p.s.storage.types.elec) for ts in p.time_steps) + 
+
+                    sum(sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
+                         for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                  )
             )
         else
             binWHL = @variable(m, binary = true)

src/constraints/generator_constraints.jl

@@ -1,60 +1,88 @@
-# REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE.
-function add_fuel_burn_constraints(m,p)
+# *********************************************************************************
+# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without modification,
+# are permitted provided that the following conditions are met:
+#
+# Redistributions of source code must retain the above copyright notice, this list
+# of conditions and the following disclaimer.
+#
+# Redistributions in binary form must reproduce the above copyright notice, this
+# list of conditions and the following disclaimer in the documentation and/or other
+# materials provided with the distribution.
+#
+# Neither the name of the copyright holder nor the names of its contributors may be
+# used to endorse or promote products derived from this software without specific
+# prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+# OF THE POSSIBILITY OF SUCH DAMAGE.
+# *********************************************************************************
+function add_fuel_burn_constraints(m,p; _n="")
 	fuel_slope_gal_per_kwhe, fuel_intercept_gal_per_hr = fuel_slope_and_intercept(
 		electric_efficiency_full_load=p.s.generator.electric_efficiency_full_load, 
 		electric_efficiency_half_load=p.s.generator.electric_efficiency_half_load,
 		fuel_higher_heating_value_kwh_per_unit=p.s.generator.fuel_higher_heating_value_kwh_per_gal
 	)
   	@constraint(m, [t in p.techs.gen, ts in p.time_steps],
-		m[:dvFuelUsage][t, ts] == (fuel_slope_gal_per_kwhe * p.s.generator.fuel_higher_heating_value_kwh_per_gal *
-		p.production_factor[t, ts] * p.hours_per_time_step * m[:dvRatedProduction][t, ts]) +
-		(fuel_intercept_gal_per_hr * p.s.generator.fuel_higher_heating_value_kwh_per_gal * p.hours_per_time_step * m[:binGenIsOnInTS][t, ts])
+		m[Symbol("dvFuelUsage"*_n)][t, ts] == (fuel_slope_gal_per_kwhe * p.s.generator.fuel_higher_heating_value_kwh_per_gal *
+		p.production_factor[t, ts] * p.hours_per_time_step * m[Symbol("dvRatedProduction"*_n)][t, ts]) +
+		(fuel_intercept_gal_per_hr * p.s.generator.fuel_higher_heating_value_kwh_per_gal * p.hours_per_time_step * m[Symbol("binGenIsOnInTS"*_n)][t, ts])
 	)
 	@constraint(m,
-		sum(m[:dvFuelUsage][t, ts] for t in p.techs.gen, ts in p.time_steps) <=
+		sum(m[Symbol("dvFuelUsage"*_n)][t, ts] for t in p.techs.gen, ts in p.time_steps) <=
 		p.s.generator.fuel_avail_gal * p.s.generator.fuel_higher_heating_value_kwh_per_gal
 	)
 end
 
 
-function add_binGenIsOnInTS_constraints(m,p)
+function add_binGenIsOnInTS_constraints(m,p; _n="")
 	# Generator must be on for nonnegative output
 	@constraint(m, [t in p.techs.gen, ts in p.time_steps],
-		m[:dvRatedProduction][t, ts] <= p.max_sizes[t] * m[:binGenIsOnInTS][t, ts]
+		m[Symbol("dvRatedProduction"*_n)][t, ts] <= p.max_sizes[t] * m[Symbol("binGenIsOnInTS"*_n)][t, ts]
 	)
 	# Note: min_turn_down_fraction is only enforced when `off_grid_flag` is true and in p.time_steps_with_grid, but not for grid outages for on-grid analyses
 	if p.s.settings.off_grid_flag 
 		@constraint(m, [t in p.techs.gen, ts in p.time_steps_without_grid],
-			p.s.generator.min_turn_down_fraction * m[:dvSize][t] - m[:dvRatedProduction][t, ts] <=
-			p.max_sizes[t] * (1 - m[:binGenIsOnInTS][t, ts])
+			p.s.generator.min_turn_down_fraction * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <=
+			p.max_sizes[t] * (1 - m[Symbol("binGenIsOnInTS"*_n)][t, ts])
 		)
 	else 
 		@constraint(m, [t in p.techs.gen, ts in p.time_steps_with_grid],
-			p.s.generator.min_turn_down_fraction * m[:dvSize][t] - m[:dvRatedProduction][t, ts] <=
-			p.max_sizes[t] * (1 - m[:binGenIsOnInTS][t, ts])
+			p.s.generator.min_turn_down_fraction * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <=
+			p.max_sizes[t] * (1 - m[Symbol("binGenIsOnInTS"*_n)][t, ts])
 		)
 	end 
 end
 
 
-function add_gen_can_run_constraints(m,p)
+function add_gen_can_run_constraints(m,p; _n="")
 	if p.s.generator.only_runs_during_grid_outage
 		for ts in p.time_steps_with_grid, t in p.techs.gen
-			fix(m[:dvRatedProduction][t, ts], 0.0, force=true)
+			fix(m[Symbol("dvRatedProduction"*_n)][t, ts], 0.0, force=true)
 		end
 	end
 
 	if !(p.s.generator.sells_energy_back_to_grid)
 		for t in p.techs.gen, u in p.export_bins_by_tech[t], ts in p.time_steps
-			fix(m[:dvProductionToGrid][t, u, ts], 0.0, force=true)
+			fix(m[Symbol("dvProductionToGrid"*_n)][t, u, ts], 0.0, force=true)
 		end
 	end
 end
 
 
-function add_gen_rated_prod_constraint(m, p)
+function add_gen_rated_prod_constraint(m, p; _n="")
 	@constraint(m, [t in p.techs.gen, ts in p.time_steps],
-		m[:dvSize][t] >= m[:dvRatedProduction][t, ts]
+		m[Symbol("dvSize"*_n)][t] >= m[Symbol("dvRatedProduction"*_n)][t, ts]
 	)
 end
 
@@ -64,17 +92,17 @@ end
 
 Add Generator operational constraints and cost expressions.
 """
-function add_gen_constraints(m, p)
-    add_fuel_burn_constraints(m,p)
-    add_binGenIsOnInTS_constraints(m,p)
-    add_gen_can_run_constraints(m,p)
-    add_gen_rated_prod_constraint(m,p)
+function add_gen_constraints(m, p; _n="")
+    add_fuel_burn_constraints(m,p, _n=_n)
+    add_binGenIsOnInTS_constraints(m,p, _n=_n)
+    add_gen_can_run_constraints(m,p, _n=_n)
+    add_gen_rated_prod_constraint(m,p, _n=_n)
 
-    m[:TotalGenPerUnitProdOMCosts] = @expression(m, p.third_party_factor * p.pwf_om *
+    m[Symbol("TotalGenPerUnitProdOMCosts"*_n)] = @expression(m, p.third_party_factor * p.pwf_om *
         sum(p.s.generator.om_cost_per_kwh * p.hours_per_time_step *
-        m[:dvRatedProduction][t, ts] for t in p.techs.gen, ts in p.time_steps)
-    )
-    m[:TotalGenFuelCosts] = @expression(m,
-        sum(p.pwf_fuel[t] * m[:dvFuelUsage][t,ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.gen, ts in p.time_steps)
+        m[Symbol("dvRatedProduction"*_n)][t, ts] for t in p.techs.gen, ts in p.time_steps)
     )
+    m[Symbol("TotalGenFuelCosts"*_n)] = @expression(m, 
+		sum(p.pwf_fuel[t] * m[Symbol("dvFuelUsage"*_n)][t,ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.gen, ts in p.time_steps) 
+	)
 end

src/constraints/storage_constraints.jl

@@ -68,7 +68,7 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
         m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) 
             + p.s.storage.attr[b].grid_charge_efficiency * m[Symbol("dvGridToStorage"*_n)][b, ts] 
-            - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency
+            - ((m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + m[Symbol("dvStorageToGrid"*_n)][b,ts])/ p.s.storage.attr[b].discharge_efficiency)
         )
 	)
 	# Constraint (4g)-2: state-of-charge for electrical storage - no grid
@@ -102,32 +102,45 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 
 	#Constraint (4k)-alt: Dispatch to and from electrical storage is no greater than power capacity
 	@constraint(m, [ts in p.time_steps_with_grid],
-        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] + 
+        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] + m[Symbol("dvStorageToGrid"*_n)][b,ts] +
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) + m[Symbol("dvGridToStorage"*_n)][b, ts]
     )
 
+	#Dispatch from electrical storage is no greater than power capacity 
+    @constraint(m, [ts in p.time_steps_without_grid],
+    m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + m[Symbol("dvStorageToGrid"*_n)][b,ts])
+
 	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity (no grid connection)
 	@constraint(m, [ts in p.time_steps_without_grid],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + 
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec)
     )
 
-    # Remove grid-to-storage as an option if option to grid charge is turned off
+    #Constraint (4m)-1: Remove grid-to-storage as an option if option to grid charge is turned off
     if !(p.s.storage.attr[b].can_grid_charge)
         for ts in p.time_steps_with_grid
             fix(m[Symbol("dvGridToStorage"*_n)][b, ts], 0.0, force=true)
         end
 	end
 
+    #Constraint (4m)-2: Force storage export to grid to zero if option to grid export is turned off
+    if !p.s.storage.attr[b].can_export_to_grid
+        for ts in p.time_steps
+            fix(m[Symbol("dvStorageToGrid"*_n)][b, ts], 0.0, force=true)
+        end
+    end
+
     if p.s.storage.attr[b].minimum_avg_soc_fraction > 0
         avg_soc = sum(m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps) /
                    (8760. / p.hours_per_time_step)
         @constraint(m, avg_soc >= p.s.storage.attr[b].minimum_avg_soc_fraction * 
             sum(m[Symbol("dvStorageEnergy"*_n)][b])
         )
     end
+
 end
 
+
 function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
 
     # Constraint (4j)-1: Reconcile state-of-charge for (hot) thermal storage

src/core/electric_utility.jl

@@ -558,7 +558,7 @@ function cambium_emissions_profile(; scenario::String,
                                     emissions_year::Int=2017
                                     )
 
-    url = "https://scenarioviewer.nrel.gov/api/get-levelized/" # Production 
+    url = "https://scenarioviewer.nlr.gov/api/get-levelized/" # Production 
     project_uuid = "82460f06-548c-4954-b2d9-b84ba92d63e2" # Cambium 2022 
 
     payload=Dict(