Merge pull request #101 from NREL-Sienna/jd/fix_modified_RTS

Jd/fix modified rts

Project.toml

@@ -1,7 +1,7 @@
 name = "PowerSystemCaseBuilder"
 uuid = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 authors = ["Sourabh Dalvi", "Jose Daniel Lara"]
-version = "1.3.1"
+version = "1.3.2"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -29,7 +29,7 @@ HDF5 = "0.17"
 InfrastructureSystems = "2"
 JSON3 = "1"
 LazyArtifacts = "1"
-PowerSystems = "4"
+PowerSystems = "^4.1.1"
 PrettyTables = "2"
 Random = "1"
 SHA = "0.7"

src/library/psi_library.jl

@@ -1418,40 +1418,47 @@ function _duplicate_system(main_sys::PSY.System, twin_sys::PSY.System, HVDC_line
     # cost perturbation must be the same for each sub-system
     rand_ix = 1
     for g in get_components(
-        x -> x.prime_mover_type in [PrimeMovers.CT, PrimeMovers.CC],
+        x -> get_fuel(x) in [ThermalFuels.NATURAL_GAS, ThermalFuels.COAL],
         ThermalStandard,
         main_sys,
     )
-        noise_values = rand(MersenneTwister(COST_PERTURBATION_NOISE_SEED), 1_000_000)
-        old_pwl_array = get_points(get_value_curve(get_variable(get_operation_cost(g))))
-        new_pwl_array = similar(old_pwl_array)
-        for (ix, (x, y)) in enumerate(old_pwl_array)
-            if ix ∈ [1, length(old_pwl_array)]
+        # This makes the twin system cheaper for the first tranche
+        # Creates an imbalance in which side is more expensive for testing
+        # purposes
+        direction = occursin("twin", PSY.get_name(g)) ? -1 : 1
+        @show PSY.get_name(g), direction
+        noise_values = rand(MersenneTwister(COST_PERTURBATION_NOISE_SEED), 10_000_000)
+        old_value_curve = get_value_curve(get_variable(get_operation_cost(g)))
+        old_slopes = get_slopes(old_value_curve)
+        new_slopes = zeros(size(old_slopes))
+        noise_val, rand_ix = iterate(noise_values, rand_ix)
+        cost_noise = round(100.0 * noise_val; digits = 2)
+        get_initial_input(old_value_curve)
+        old_y =
+            get_initial_input(old_value_curve) /
+            (get_active_power_limits(g).min * get_base_power(g))
+        new_first_input =
+            (old_y + direction * cost_noise) * get_active_power_limits(g).min *
+            get_base_power(g)
+        new_slopes[1] = old_slopes[1] + direction * cost_noise
+        @assert new_slopes[1] > 0.0
+        for ix in 2:length(old_slopes)
+            while new_slopes[ix - 1] > new_slopes[ix]
                 noise_val, rand_ix = iterate(noise_values, rand_ix)
-                cost_noise = 50.0 * noise_val
-                new_pwl_array[ix] = (x = x, y = (y + cost_noise))
-            else
-                try_again = true
-                while try_again
-                    noise_val, rand_ix = iterate(noise_values, rand_ix)
-                    cost_noise = 50.0 * noise_val
-                    noise_val, rand_ix = iterate(noise_values, rand_ix)
-                    power_noise = 0.01 * noise_val
-                    slope_previous =
-                        ((y + cost_noise) - old_pwl_array[ix - 1].y) /
-                        ((x - power_noise) - old_pwl_array[ix - 1].x)
-                    slope_next =
-                        (-(y + cost_noise) + old_pwl_array[ix + 1].y) /
-                        (-(x - power_noise) + old_pwl_array[ix + 1].x)
-                    new_pwl_array[ix] = (x = (x - power_noise), y = (y + cost_noise))
-                    try_again = slope_previous > slope_next
-                    if rand_ix == length(noise_values)
-                        break
-                    end
-                end
+                cost_noise = round(100.0 * noise_val; digits = 2)
+                new_slopes[ix] = old_slopes[ix] + cost_noise
             end
         end
-        set_variable!(get_operation_cost(g), CostCurve(PiecewisePointCurve(new_pwl_array)))
+        @assert old_slopes != new_slopes
+        set_variable!(
+            get_operation_cost(g),
+            CostCurve(
+                PiecewiseIncrementalCurve(
+                    nothing,
+                    new_first_input,
+                    get_x_coords(old_value_curve),
+                    new_slopes,
+                )))
     end
 
     # set service participation