Merge branch 'develop' into rate-bug-fix

CHANGELOG.md

@@ -26,6 +26,16 @@ Classify the change according to the following categories:
 ## Develop 
 ### Fixed
 - Don't double add adjustments to urdb rates with non-standard units 
+
+## Develop - 2023-06-21
+### Changed
+- Consolidated PVWatts API calls to 1 call (previously 3 separate calls existed). API call occurs in `src/core/utils.jl/call_pvwatts_api()`. This function is called for PV in `src/core/production_factor.jl/get_production_factor(PV)` and for GHP in `src/core/scenario.jl`. If GHP and PV are evaluated together, the GHP PVWatts call for ambient temperature is also used to assign the pv.production_factor_series in Scenario.jl so that the PVWatts API does not get called again downstream in `get_production_factor(PV)`.  
+- In `src/core/utils.jl/call_pvwatts_api()`, updated NSRDB bounds used in PVWatts query (now includes southern New Zealand)
+- Updated PV Watts version from v6 to v8. PVWatts V8 updates the weather data to 2020 TMY data from the NREL NSRDB for locations covered by the database. (The NSRDB weather data used in PVWatts V6 is from around 2015.) See other differences at https://developer.nrel.gov/docs/solar/pvwatts/.
+- Made PV struct mutable: This allows for assigning pv.production_factor_series when calling PVWatts for GHP, to avoid a extra PVWatts calls later.
+### Fixed
+- Issue with using a leap year with a URDB rate - the URDB rate was creating energy_rate of length 8784 instead of intended 8760
+
 ## v0.32.3
 ### Fixed
 - Calculate **num_battery_bins** default in `backup_reliability.jl` based on battery duration to prevent significant discretization error (and add test)

src/core/production_factor.jl

@@ -35,44 +35,13 @@ function get_production_factor(pv::PV, latitude::Real, longitude::Real; timefram
         return pv.production_factor_series
     end
 
-    # Check if site is beyond the bounds of the NRSDB dataset. If so, use the international dataset.
-    dataset = "nsrdb"
-    if longitude < -179.5 || longitude > -21.0 || latitude < -21.5 || latitude > 60.0
-        if longitude < 81.5 || longitude > 179.5 || latitude < -43.8 || latitude > 60.0 
-            if longitude < 67.0 || longitude > 81.5 || latitude < -43.8 || latitude > 38.0
-                dataset = "intl"
-            end
-        end
-    end
+    watts, ambient_temp_celcius = call_pvwatts_api(latitude, longitude; tilt=pv.tilt, azimuth=pv.azimuth, module_type=pv.module_type, 
+        array_type=pv.array_type, losses=round(pv.losses*100, digits=3), dc_ac_ratio=pv.dc_ac_ratio,
+        gcr=pv.gcr, inv_eff=pv.inv_eff*100, timeframe=timeframe, radius=pv.radius,
+        time_steps_per_hour=time_steps_per_hour)
 
-    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
-        "&lat=", latitude , "&lon=", longitude, "&tilt=", pv.tilt,
-        "&system_capacity=1", "&azimuth=", pv.azimuth, "&module_type=", pv.module_type,
-        "&array_type=", pv.array_type, "&losses=", round(pv.losses*100, digits=3), "&dc_ac_ratio=", pv.dc_ac_ratio,
-        "&gcr=", pv.gcr, "&inv_eff=", pv.inv_eff*100, "&timeframe=", timeframe, "&dataset=", dataset,
-        "&radius=", pv.radius
-    )
+    return watts
 
-    try
-        @info "Querying PVWatts for production_factor with " pv.name
-        r = HTTP.get(url, keepalive=true, readtimeout=10)
-        @info "Response received from PVWatts"
-        response = JSON.parse(String(r.body))
-        if r.status != 200
-            throw(@error("Bad response from PVWatts: $(response["errors"])"))
-        end
-        @info "PVWatts success."
-        watts = collect(get(response["outputs"], "ac", []) / 1000)  # scale to 1 kW system (* 1 kW / 1000 W)
-        if length(watts) != 8760
-            throw(@error("PVWatts did not return a valid production factor. Got $watts"))
-        end
-        if time_steps_per_hour > 1
-            watts = repeat(watts, inner=time_steps_per_hour)
-        end
-        return watts
-    catch e
-        throw(@error("Error occurred when calling PVWatts: $e"))
-    end
 end
 
 

src/core/pv.jl

@@ -33,12 +33,12 @@
     array_type::Int=1, # PV Watts array type (0: Ground Mount Fixed (Open Rack); 1: Rooftop, Fixed; 2: Ground Mount 1-Axis Tracking; 3 : 1-Axis Backtracking; 4: Ground Mount, 2-Axis Tracking)
     tilt::Real= array_type == 1 ? 10 : abs(latitude), # tilt = 10 deg for rooftop systems, abs(lat) for ground-mount
     module_type::Int=0, # PV module type (0: Standard; 1: Premium; 2: Thin Film)
-    losses::Real=0.14,
+    losses::Real=0.14, # System losses
     azimuth::Real = latitude≥0 ? 180 : 0, # set azimuth to zero for southern hemisphere
-    gcr::Real=0.4,
-    radius::Int=0,
-    name::String="PV",
-    location::String="both",
+    gcr::Real=0.4,  # Ground coverage ratio
+    radius::Int=0, # Radius, in miles, to use when searching for the closest climate data station. Use zero to use the closest station regardless of the distance
+    name::String="PV", # for use with multiple pvs 
+    location::String="both", # one of ["roof", "ground", "both"]
     existing_kw::Real=0,
     min_kw::Real=0,
     max_kw::Real=1.0e9,
@@ -82,7 +82,7 @@
     If `azimuth` is not provided, then it is set to 180 if the site is in the northern hemisphere and 0 if in the southern hemisphere.
 
 """
-struct PV <: AbstractTech
+mutable struct PV <: AbstractTech
     tilt
     array_type
     module_type
@@ -151,7 +151,7 @@ struct PV <: AbstractTech
         acres_per_kw::Real=6e-3,
         inv_eff::Real=0.96,
         dc_ac_ratio::Real=1.2,
-        production_factor_series::Union{Nothing, Array{Real,1}} = nothing,
+        production_factor_series::Union{Nothing, Array{<:Real,1}} = nothing,
         federal_itc_fraction::Real = 0.3,
         federal_rebate_per_kw::Real = 0.0,
         state_ibi_fraction::Real = 0.0,

src/core/scenario.jl

@@ -454,37 +454,27 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
             number_of_ghpghx = length(d["GHP"]["ghpghx_inputs"])
         end
         # Call PVWatts for hourly dry-bulb outdoor air temperature
-        ambient_temperature_f = []
+        ambient_temp_degF = []
         if !haskey(d["GHP"]["ghpghx_inputs"][1], "ambient_temperature_f") || isempty(d["GHP"]["ghpghx_inputs"][1]["ambient_temperature_f"])
-            url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
-                    "&lat=", d["Site"]["latitude"] , "&lon=", d["Site"]["longitude"], "&tilt=", d["Site"]["latitude"],
-                    "&system_capacity=1", "&azimuth=", 180, "&module_type=", 0,
-                    "&array_type=", 0, "&losses=", 0.14, "&dc_ac_ratio=", 1.1,
-                    "&gcr=", 0.4, "&inv_eff=", 99, "&timeframe=", "hourly", "&dataset=nsrdb",
-                    "&radius=", 100)
-            try
-                @info "Querying PVWatts for ambient temperature"
-                r = HTTP.get(url)
-                response = JSON.parse(String(r.body))
-                if r.status != 200
-                    throw(@error("Bad response from PVWatts: $(response["errors"])"))
+            # If PV is evaluated and we need to call PVWatts for ambient temperature, assign PV production factor here too with the same call
+            # By assigning pv.production_factor_series here, it will skip the PVWatts call in get_production_factor(PV) call from reopt_input.jl
+            if !isempty(pvs)
+                for pv in pvs
+                    pv.production_factor_series, ambient_temp_celcius = call_pvwatts_api(site.latitude, site.longitude; tilt=pv.tilt, azimuth=pv.azimuth, module_type=pv.module_type, 
+                        array_type=pv.array_type, losses=round(pv.losses*100, digits=3), dc_ac_ratio=pv.dc_ac_ratio,
+                        gcr=pv.gcr, inv_eff=pv.inv_eff*100, timeframe="hourly", radius=pv.radius, time_steps_per_hour=settings.time_steps_per_hour)
                 end
-                @info "PVWatts success."
-                temp_c = get(response["outputs"], "tamb", [])
-                if length(temp_c) != 8760 || isempty(temp_c)
-                    throw(@error("PVWatts did not return a valid temperature profile. Got $temp_c"))
-                end
-                ambient_temperature_f = temp_c * 1.8 .+ 32.0
-            catch e
-                throw(@error("Error occurred when calling PVWatts: $e"))
+            else
+                pv_prodfactor, ambient_temp_celcius = call_pvwatts_api(site.latitude, site.longitude; time_steps_per_hour=settings.time_steps_per_hour)    
             end
+            ambient_temp_degF = ambient_temp_celcius * 1.8 .+ 32.0
         else
-            ambient_temperature_f = d["GHP"]["ghpghx_inputs"][1]["ambient_temperature_f"]
+            ambient_temp_degF = d["GHP"]["ghpghx_inputs"][1]["ambient_temperature_f"]
         end
 
         for i in 1:number_of_ghpghx
             ghpghx_inputs = d["GHP"]["ghpghx_inputs"][i]
-            d["GHP"]["ghpghx_inputs"][i]["ambient_temperature_f"] = ambient_temperature_f
+            d["GHP"]["ghpghx_inputs"][i]["ambient_temperature_f"] = ambient_temp_degF
             # Only SpaceHeating portion of Heating Load gets served by GHP, unless allowed by can_serve_dhw
             if get(ghpghx_inputs, "heating_thermal_load_mmbtu_per_hr", []) in [nothing, []]
                 if haskey(d["GHP"], "can_serve_dhw") && d["GHP"]["can_serve_dhw"]

src/core/site.jl

@@ -36,7 +36,7 @@ Inputs related to the physical location:
     longitude::Real, 
     land_acres::Union{Real, Nothing} = nothing, # acres of land available for PV panels and/or Wind turbines. Constraint applied separately to PV and Wind, meaning the two technologies are assumed to be able to be co-located.
     roof_squarefeet::Union{Real, Nothing} = nothing,
-    min_resil_time_steps::Int=0,
+    min_resil_time_steps::Int=0, # The minimum number consecutive timesteps that load must be fully met once an outage begins. Only applies to multiple outage modeling using inputs outage_start_time_steps and outage_durations.
     mg_tech_sizes_equal_grid_sizes::Bool = true,
     node::Int = 1,
     CO2_emissions_reduction_min_fraction::Union{Float64, Nothing} = nothing,

src/core/urdb.jl

@@ -270,6 +270,9 @@ function parse_urdb_energy_costs(d::Dict, year::Int; time_steps_per_hour=1, bigM
 
         for month in range(1, stop=12)
             n_days = daysinmonth(Date(string(year) * "-" * string(month)))
+            if month == 2 && isleapyear(year)
+                n_days -= 1
+            end
 
             for day in range(1, stop=n_days)
 

src/core/utils.jl

@@ -377,54 +377,57 @@ function generate_year_profile_hourly(year::Int64, consecutive_periods::Abstract
 end
 
 
-function get_ambient_temperature(latitude::Real, longitude::Real; timeframe="hourly")
-    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
-        "&lat=", latitude , "&lon=", longitude, "&tilt=", latitude,
-        "&system_capacity=1", "&azimuth=", 180, "&module_type=", 0,
-        "&array_type=", 0, "&losses=", 14,
-        "&timeframe=", timeframe, "&dataset=nsrdb"
-    )
-
-    try
-        @info "Querying PVWatts for ambient temperature... "
-        r = HTTP.get(url)
-        response = JSON.parse(String(r.body))
-        if r.status != 200
-            throw(@error("Bad response from PVWatts: $(response["errors"])"))
-        end
-        @info "PVWatts success."
-        tamb = collect(get(response["outputs"], "tamb", []))  # Celcius
-        if length(tamb) != 8760
-            throw(@error("PVWatts did not return a valid temperature. Got $tamb"))
+"""
+    call_pvwatts_api(latitude::Real, longitude::Real; tilt=latitude, azimuth=180, module_type=0, array_type=1, 
+        losses=14, dc_ac_ratio=1.2, gcr=0.4, inv_eff=96, timeframe="hourly", radius=0, time_steps_per_hour=1)
+This calls the PVWatts API and returns both:
+ - PV production factor
+ - Ambient outdoor air dry bulb temperature profile [Celcius] 
+"""
+function call_pvwatts_api(latitude::Real, longitude::Real; tilt=latitude, azimuth=180, module_type=0, array_type=1, 
+    losses=14, dc_ac_ratio=1.2, gcr=0.4, inv_eff=96, timeframe="hourly", radius=0, time_steps_per_hour=1)
+    # Check if site is beyond the bounds of the NRSDB TMY dataset. If so, use the international dataset.
+    dataset = "nsrdb"
+    if longitude < -179.5 || longitude > -21.0 || latitude < -21.5 || latitude > 60.0
+        if longitude < 81.5 || longitude > 179.5 || latitude < -60.0 || latitude > 60.0 
+            if longitude < 67.0 || latitude < -40.0 || latitude > 38.0
+                dataset = "intl"
+            end
         end
-        return tamb
-    catch e
-        throw(@error("Error occurred when calling PVWatts: $e"))
     end
-end
-
-
-function get_pvwatts_prodfactor(latitude::Real, longitude::Real; timeframe="hourly")
-    url = string("https://developer.nrel.gov/api/pvwatts/v6.json", "?api_key=", nrel_developer_key,
-        "&lat=", latitude , "&lon=", longitude, "&tilt=", latitude,
-        "&system_capacity=1", "&azimuth=", 180, "&module_type=", 0,
-        "&array_type=", 0, "&losses=", 14,
-        "&timeframe=", timeframe, "&dataset=nsrdb"
-    )
+    url = string("https://developer.nrel.gov/api/pvwatts/v8.json", "?api_key=", nrel_developer_key,
+        "&lat=", latitude , "&lon=", longitude, "&tilt=", tilt,
+        "&system_capacity=1", "&azimuth=", azimuth, "&module_type=", module_type,
+        "&array_type=", array_type, "&losses=", losses, "&dc_ac_ratio=", dc_ac_ratio,
+        "&gcr=", gcr, "&inv_eff=", inv_eff, "&timeframe=", timeframe, "&dataset=", dataset,
+        "&radius=", radius
+        )
 
     try
-        @info "Querying PVWatts for production factor of 1 kW system with tilt set to latitude... "
-        r = HTTP.get(url)
+        @info "Querying PVWatts for production factor and ambient air temperature... "
+        r = HTTP.get(url, keepalive=true, readtimeout=10)
         response = JSON.parse(String(r.body))
         if r.status != 200
             throw(@error("Bad response from PVWatts: $(response["errors"])"))
         end
         @info "PVWatts success."
+        # Get both possible data of interest
         watts = collect(get(response["outputs"], "ac", []) / 1000)  # scale to 1 kW system (* 1 kW / 1000 W)
+        tamb_celcius = collect(get(response["outputs"], "tamb", []))  # Celcius
+        # Validate outputs
         if length(watts) != 8760
             throw(@error("PVWatts did not return a valid prodfactor. Got $watts"))
         end
-        return watts
+        # Validate tamb_celcius
+        if length(tamb_celcius) != 8760
+            throw(@error("PVWatts did not return a valid temperature. Got $tamb_celcius"))
+        end 
+        # Upsample or downsample based on model time_steps_per_hour
+        if time_steps_per_hour > 1
+            watts = repeat(watts, inner=time_steps_per_hour)
+            tamb_celcius = repeat(tamb_celcius, inner=time_steps_per_hour)
+        end
+        return watts, tamb_celcius
     catch e
         throw(@error("Error occurred when calling PVWatts: $e"))
     end

test/runtests.jl

@@ -77,8 +77,8 @@ else  # run HiGHS tests
         inputs = REoptInputs(s)
         results = run_reopt(model, inputs)
 
-        @test results["PV"]["size_kw"] ≈ 70.3084 atol=0.01
-        @test results["Financial"]["lcc"] ≈ 430747.0 rtol=1e-5 # with levelization_factor hack the LCC is within 5e-5 of REopt API LCC
+        @test results["PV"]["size_kw"] ≈ 68.9323 atol=0.01
+        @test results["Financial"]["lcc"] ≈ 432672.0 rtol=1e-5 # with levelization_factor hack the LCC is within 5e-5 of REopt API LCC
         @test all(x == 0.0 for x in results["PV"]["electric_to_load_series_kw"][1:744])
     end
 
@@ -98,17 +98,17 @@ else  # run HiGHS tests
         r = run_reopt(model, "./scenarios/pv_storage.json")
 
         @test r["PV"]["size_kw"] ≈ 216.6667 atol=0.01
-        @test r["Financial"]["lcc"] ≈ 1.240037e7 rtol=1e-5
-        @test r["ElectricStorage"]["size_kw"] ≈ 55.9 atol=0.1
-        @test r["ElectricStorage"]["size_kwh"] ≈ 78.9 atol=0.1
+        @test r["Financial"]["lcc"] ≈ 1.239151e7 rtol=1e-5
+        @test r["ElectricStorage"]["size_kw"] ≈ 49.0 atol=0.1
+        @test r["ElectricStorage"]["size_kwh"] ≈ 83.3 atol=0.1
     end
 
     @testset "Outage with Generator" begin
         model = Model(optimizer_with_attributes(HiGHS.Optimizer, 
             "output_flag" => false, "log_to_console" => false)
         )
         results = run_reopt(model, "./scenarios/generator.json")
-        @test results["Generator"]["size_kw"] ≈ 8.13 atol=0.01
+        @test results["Generator"]["size_kw"] ≈ 9.53 atol=0.01
         @test (sum(results["Generator"]["electric_to_load_series_kw"][i] for i in 1:9) + 
             sum(results["Generator"]["electric_to_load_series_kw"][i] for i in 13:8760)) == 0
         p = REoptInputs("./scenarios/generator.json")
@@ -464,7 +464,7 @@ else  # run HiGHS tests
 
         @test reliability_results["unlimited_fuel_cumulative_survival_final_time_step"][1] ≈ 0.802997 atol=0.0001
         @test reliability_results["cumulative_survival_final_time_step"][1] ≈ 0.802997 atol=0.0001
-        @test reliability_results["mean_cumulative_survival_final_time_step"] ≈ 0.817088 atol=0.0001
+        @test reliability_results["mean_cumulative_survival_final_time_step"] ≈ 0.817586 atol=0.0001
     end                            
 
     # removed Wind test for two reasons

test/test_with_cplex.jl

@@ -59,9 +59,9 @@ end
     results = run_reopt(model, "./scenarios/pv_storage.json")
 
     @test results["PV"]["size_kw"] ≈ 217 atol=1
-    @test results["Financial"]["lcc"] ≈ 1.240037e7 rtol=1e-5
-    @test results["ElectricStorage"]["size_kw"] ≈ 56 atol=1
-    @test results["ElectricStorage"]["size_kwh"] ≈ 79 atol=1
+    @test results["Financial"]["lcc"] ≈ 1.239151e7 rtol=1e-5
+    @test results["ElectricStorage"]["size_kw"] ≈ 49 atol=1
+    @test results["ElectricStorage"]["size_kwh"] ≈ 83 atol=1
 end
 
 
@@ -74,7 +74,7 @@ end
     @test value(m[:binMGTechUsed]["Generator"]) == 1
     @test value(m[:binMGTechUsed]["PV"]) == 1
     @test value(m[:binMGStorageUsed]) == 1
-    @test results["Financial"]["lcc"] ≈ 7.19753998668e7 atol=5e4
+    @test results["Financial"]["lcc"] ≈ 6.82164056207e7 atol=5e4
 
     #=
     Scenario with $0/kWh value_of_lost_load_per_kwh, 12x169 hour outages, 1kW load/hour, and min_resil_time_steps = 168
@@ -98,7 +98,7 @@ end
         REoptInputs("./scenarios/monthly_rate.json"),
     ];
     results = run_reopt(m, ps)
-    @test results[3]["Financial"]["lcc"] + results[10]["Financial"]["lcc"] ≈ 1.23887e7 + 437169.0 rtol=1e-5
+    @test results[3]["Financial"]["lcc"] + results[10]["Financial"]["lcc"] ≈ 1.2830591384e7 rtol=1e-5
 end