Additional constraints to prevent violation of state of charge limits…

… in non-representative periods (#781)
GenXProject · Nov 18, 2024 · 6fabc96 · 6fabc96
1 parent 65e25d7
commit 6fabc96
Show file tree

Hide file tree

Showing 5 changed files with 183 additions and 2 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fusion plant optional features for thermal plants (#743).
 - Support for reusing the same Gurobi environment for multiple solves when 
 number of concurrent Gurobi uses is limited (#783).
+- Additional long-duration storage constraints to bound state of charge in 
+non-representative periods (#781).
 
 ### Changed
 - The `charge.csv` and `storage.csv` files now include only resources with 

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "GenX"
 uuid = "5d317b1e-30ec-4ed6-a8ce-8d2d88d7cfac"
 authors = ["Bonaldo, Luca", "Chakrabarti, Sambuddha", "Cheng, Fangwei", "Ding, Yifu", "Jenkins, Jesse D.", "Luo, Qian", "Macdonald, Ruaridh", "Mallapragada, Dharik", "Manocha, Aneesha", "Mantegna, Gabe ", "Morris, Jack", "Patankar, Neha", "Pecci, Filippo", "Schwartz, Aaron", "Schwartz, Jacob", "Schivley, Greg", "Sepulveda, Nestor", "Xu, Qingyu", "Zhou, Justin"]
-version = "0.4.1-dev.10"
+version = "0.4.1-dev.11"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

diff --git a/src/model/resources/hydro/hydro_inter_period_linkage.jl b/src/model/resources/hydro/hydro_inter_period_linkage.jl
@@ -42,6 +42,43 @@ Finally, the next constraint enforces that the initial storage level for each in
 \quad \forall n \in \mathcal{N}, o \in \mathcal{O}^{LDES}
 \end{aligned}
 ```
+
+**Bound storage inventory in non-representative periods**
+We need additional variables and constraints to ensure that the storage content is within zero and the installed energy capacity in non-representative periods. We introduce
+the variables $\Delta Q^{max,pos}_{o,z,m}$ and $\Delta Q^{max,neg}_{o,z,m}$ that represent the maximum positive and negative storage content variations within the representative
+period $m$, respectively, extracted as:
+
+```math
+\begin{aligned}
+& \Delta Q^{max,pos}_{o,z,m} \geq \Gamma_{o,z,(m-1)\times \tau^{period}+t } - \Gamma_{o,z,(m-1)\times \tau^{period}+1 } \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}, t \in \mathcal{T}
+& \end{aligned}
+```
+
+```math
+\begin{aligned}
+& \Delta Q^{max,neg}_{o,z,m} \leq \Gamma_{o,z,(m-1)\times \tau^{period}+t } - \Gamma_{o,z,(m-1)\times \tau^{period}+1 } \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}, t \in \mathcal{T}
+& \end{aligned}
+```
+
+For every input period $n \in \mathcal{N}$, the maximum storage is computed and constrained to be less than or equal to the installed energy capacity as:
+
+```math
+\begin{aligned}
+&  Q_{o,z,n} \times \left(1-\eta_{o,z}^{loss}\right) - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,(m-1)\times \tau^{period}+1} + \eta_{o,z}^{charge}\Pi_{o,z,(m-1)\times \tau^{period}+1} + \Delta Q^{max,pos}_{o,z,f(n)} \leq \Delta^{total, energy}_{o,z} \\
+& \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in \mathcal{N}
+& \end{aligned}
+```
+
+Similarly, the minimum storage content is imposed to be positive in every period of the time horizon:
+
+```math
+\begin{aligned}
+&  Q_{o,z,n} \times \left(1-\eta_{o,z}^{loss}\right) - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,(m-1)\times \tau^{period}+1} + \eta_{o,z}^{charge}\Pi_{o,z,(m-1)\times \tau^{period}+1} + \Delta Q^{max,neg}_{o,z,f(n)} \geq 0 \\
+& \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in \mathcal{N}
+& \end{aligned}
+```
+
+Additional details on this approach are available in [Parolin et al., 2024](https://doi.org/10.48550/arXiv.2409.19079).
 """
 function hydro_inter_period_linkage!(EP::Model, inputs::Dict)
     println("Long Duration Storage Module for Hydro Reservoir")
@@ -71,6 +108,12 @@ function hydro_inter_period_linkage!(EP::Model, inputs::Dict)
     # Build up inventory can be positive or negative
     @variable(EP, vdSOC_HYDRO[y in STOR_HYDRO_LONG_DURATION, w = 1:REP_PERIOD])
 
+    # Maximum positive storage inventory change within subperiod
+	@variable(EP, vdSOC_maxPos_HYDRO[y in STOR_HYDRO_LONG_DURATION, w=1:REP_PERIOD] >= 0)
+
+	# Maximum negative storage inventory change within subperiod
+	@variable(EP, vdSOC_maxNeg_HYDRO[y in STOR_HYDRO_LONG_DURATION, w=1:REP_PERIOD] <= 0)	
+
     ### Constraints ###
 
     # Links last time step with first time step, ensuring position in hour 1 is within eligible change from final hour position
@@ -111,4 +154,27 @@ function hydro_inter_period_linkage!(EP::Model, inputs::Dict)
             r in REP_PERIODS_INDEX],
         vSOC_HYDROw[y,r]==EP[:vS_HYDRO][y, hours_per_subperiod * dfPeriodMap[r, :Rep_Period_Index]] -
                 vdSOC_HYDRO[y, dfPeriodMap[r, :Rep_Period_Index]])
+
+    # Extract maximum storage level variation (positive) within subperiod
+    @constraint(EP, cMaxSoCVarPos_H[y in STOR_HYDRO_LONG_DURATION, w=1:REP_PERIOD, t=2:hours_per_subperiod],
+                    vdSOC_maxPos_HYDRO[y,w] >= EP[:vS_HYDRO][y,hours_per_subperiod*(w-1)+t] - EP[:vS_HYDRO][y,hours_per_subperiod*(w-1)+1])
+
+    # Extract maximum storage level variation (negative) within subperiod
+    @constraint(EP, cMaxSoCVarNeg_H[y in STOR_HYDRO_LONG_DURATION, w=1:REP_PERIOD, t=2:hours_per_subperiod],
+                    vdSOC_maxNeg_HYDRO[y,w] <= EP[:vS_HYDRO][y,hours_per_subperiod*(w-1)+t] - EP[:vS_HYDRO][y,hours_per_subperiod*(w-1)+1])
+
+    # Max storage content within each modeled period cannot exceed installed energy capacity
+    @constraint(EP, cSoCLongDurationStorageMaxInt_H[y in STOR_HYDRO_LONG_DURATION, r in MODELED_PERIODS_INDEX],
+        vSOC_HYDROw[y,r]-(1/efficiency_down(gen[y])*EP[:vP][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+        -EP[:vSPILL][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1]
+        +inputs["pP_Max"][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1]*EP[:eTotalCap][y]
+        +vdSOC_maxPos_HYDRO[y,dfPeriodMap[r,:Rep_Period_Index]] <= hydro_energy_to_power_ratio(gen[y])*EP[:eTotalCap][y])
+
+    # Min storage content within each modeled period cannot be negative
+    @constraint(EP, cSoCLongDurationStorageMinInt_H[y in STOR_HYDRO_LONG_DURATION, r in MODELED_PERIODS_INDEX],
+        vSOC_HYDROw[y,r]-(1/efficiency_down(gen[y])*EP[:vP][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+        -EP[:vSPILL][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1]
+        +inputs["pP_Max"][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1]*EP[:eTotalCap][y]
+        +vdSOC_maxPos_HYDRO[y,dfPeriodMap[r,:Rep_Period_Index]] >= 0)     
+
 end
diff --git a/src/model/resources/storage/long_duration_storage.jl b/src/model/resources/storage/long_duration_storage.jl
@@ -55,7 +55,44 @@ If the capacity reserve margin constraint is enabled, a similar set of constrain
 & \frac{1}{\eta_{o,z}^{discharge}}\Theta^{CRM}_{o,z,(m-1)\times \tau^{period}+1} - \eta_{o,z}^{charge}\Pi^{CRM}_{o,z,(m-1)\times \tau^{period}+1} \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}
 \end{aligned}
 ```
-All other constraints are identical to those used to track the actual state of charge, except with the new variables $Q^{CRM}_{o,z,n}$ and $\Delta Q^{CRM}_{o,z,n}$ used in place of $Q_{o,z,n}$ and $\Delta Q_{o,z,n}$, respectively.
+All other constraints are identical to those used to track the actual state of charge, except with the new variables $Q^{CRM}_{o,z,n}$ and $\Delta Q^{CRM}_{o,z,n}$ used in place of $Q_{o,z,n}$ and $\Delta Q_{o,z,n}$, respectively. \
+
+**Bound storage inventory in non-representative periods**
+We need additional variables and constraints to ensure that the storage content is within zero and the installed energy capacity in non-representative periods. We introduce
+the variables $\Delta Q^{max,pos}_{o,z,m}$ and $\Delta Q^{max,neg}_{o,z,m}$ that represent the maximum positive and negative storage content variations within the representative
+period $m$, respectively, extracted as:
+
+```math
+\begin{aligned}
+& \Delta Q^{max,pos}_{o,z,m} \geq \Gamma_{o,z,(m-1)\times \tau^{period}+t } - \Gamma_{o,z,(m-1)\times \tau^{period}+1 } \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}, t \in \mathcal{T}
+& \end{aligned}
+```
+
+```math
+\begin{aligned}
+& \Delta Q^{max,neg}_{o,z,m} \leq \Gamma_{o,z,(m-1)\times \tau^{period}+t } - \Gamma_{o,z,(m-1)\times \tau^{period}+1 } \quad \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, m \in \mathcal{M}, t \in \mathcal{T}
+& \end{aligned}
+```
+
+For every input period $n \in \mathcal{N}$, the maximum storage is computed and constrained to be less than or equal to the installed energy capacity as:
+
+```math
+\begin{aligned}
+&  Q_{o,z,n} \times \left(1-\eta_{o,z}^{loss}\right) - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,(m-1)\times \tau^{period}+1} + \eta_{o,z}^{charge}\Pi_{o,z,(m-1)\times \tau^{period}+1} + \Delta Q^{max,pos}_{o,z,f(n)} \leq \Delta^{total, energy}_{o,z} \\
+& \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in \mathcal{N}
+& \end{aligned}
+```
+
+Similarly, the minimum storage content is imposed to be positive in every period of the time horizon:
+
+```math
+\begin{aligned}
+&  Q_{o,z,n} \times \left(1-\eta_{o,z}^{loss}\right) - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,(m-1)\times \tau^{period}+1} + \eta_{o,z}^{charge}\Pi_{o,z,(m-1)\times \tau^{period}+1} + \Delta Q^{max,neg}_{o,z,f(n)} \geq 0 \\
+& \forall o \in \mathcal{O}^{LDES}, z \in \mathcal{Z}, n \in \mathcal{N}
+& \end{aligned}
+```
+
+Additional details on this approach are available in [Parolin et al., 2024](https://doi.org/10.48550/arXiv.2409.19079).
 """
 function long_duration_storage!(EP::Model, inputs::Dict, setup::Dict)
     println("Long Duration Storage Module")
@@ -96,6 +133,12 @@ function long_duration_storage!(EP::Model, inputs::Dict, setup::Dict)
         @variable(EP, vCAPRES_dsoc[y in STOR_LONG_DURATION, w = 1:REP_PERIOD])
     end
 
+    # Maximum positive storage inventory change within subperiod
+	@variable(EP, vdSOC_maxPos[y in STOR_LONG_DURATION, w=1:REP_PERIOD] >= 0)
+
+	# Maximum negative storage inventory change within subperiod
+	@variable(EP, vdSOC_maxNeg[y in STOR_LONG_DURATION, w=1:REP_PERIOD] <= 0)
+
     ### Constraints ###
 
     # Links last time step with first time step, ensuring position in hour 1 is within eligible change from final hour position
@@ -181,4 +224,24 @@ function long_duration_storage!(EP::Model, inputs::Dict, setup::Dict)
                 r in MODELED_PERIODS_INDEX],
             vSOCw[y, r]>=vCAPRES_socw[y, r])
     end
+
+    # Extract maximum storage level variation (positive) within subperiod
+	@constraint(EP, cMaxSoCVarPos[y in STOR_LONG_DURATION, w=1:REP_PERIOD, t=2:hours_per_subperiod],
+                    vdSOC_maxPos[y,w] >= EP[:vS][y,hours_per_subperiod*(w-1)+t] - EP[:vS][y,hours_per_subperiod*(w-1)+1])
+
+    # Extract maximum storage level variation (negative) within subperiod
+    @constraint(EP, cMaxSoCVarNeg[y in STOR_LONG_DURATION, w=1:REP_PERIOD, t=2:hours_per_subperiod],
+                    vdSOC_maxNeg[y,w] <= EP[:vS][y,hours_per_subperiod*(w-1)+t] - EP[:vS][y,hours_per_subperiod*(w-1)+1])
+
+    # Max storage content within each modeled period cannot exceed installed energy capacity
+    @constraint(EP, cSoCLongDurationStorageMaxInt[y in STOR_LONG_DURATION, r in MODELED_PERIODS_INDEX],
+            (1-self_discharge(gen[y]))*vSOCw[y,r]-(1/efficiency_down(gen[y])*EP[:vP][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+            +(efficiency_up(gen[y])*EP[:vCHARGE][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+            +vdSOC_maxPos[y,dfPeriodMap[r,:Rep_Period_Index]] <= EP[:eTotalCapEnergy][y])
+
+    # Min storage content within each modeled period cannot be negative
+    @constraint(EP, cSoCLongDurationStorageMinInt[y in STOR_LONG_DURATION, r in MODELED_PERIODS_INDEX],
+            (1-self_discharge(gen[y]))*vSOCw[y,r]-(1/efficiency_down(gen[y])*EP[:vP][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+            +(efficiency_up(gen[y])*EP[:vCHARGE][y,hours_per_subperiod*(dfPeriodMap[r,:Rep_Period_Index]-1)+1])
+            +vdSOC_maxNeg[y,dfPeriodMap[r,:Rep_Period_Index]] >= 0)
 end
diff --git a/src/write_outputs/long_duration_storage/write_opwrap_lds_stor_init.jl b/src/write_outputs/long_duration_storage/write_opwrap_lds_stor_init.jl
@@ -32,4 +32,54 @@ function write_opwrap_lds_stor_init(path::AbstractString,
     CSV.write(joinpath(path, "StorageInit.csv"),
         dftranspose(dfStorageInit, false),
         header = false)
+
+    # Write storage evolution over full time horizon
+    hours_per_subperiod = inputs["hours_per_subperiod"];
+    t_interior = 2:hours_per_subperiod
+    T_hor = hours_per_subperiod*NPeriods # total number of time steps in time horizon
+    SOC_t = zeros(G, T_hor)
+    stor_long_duration = inputs["STOR_LONG_DURATION"]
+    stor_hydro_long_duration = inputs["STOR_HYDRO_LONG_DURATION"]
+    period_map = inputs["Period_Map"].Rep_Period_Index
+    pP_max = inputs["pP_Max"]
+    e_total_cap = value.(EP[:eTotalCap])
+    v_charge = value.(EP[:vCHARGE])
+    v_P = value.(EP[:vP])
+    if setup["ParameterScale"] == 1
+        v_charge *= ModelScalingFactor
+        v_P *= ModelScalingFactor
+    end
+    if !isempty(stor_hydro_long_duration)
+        v_spill = value.(EP[:vSPILL])
+    end
+    for r in 1:NPeriods
+        w = period_map[r]
+        t_r = hours_per_subperiod * (r - 1) + 1
+        t_start_w = hours_per_subperiod * (w - 1) + 1
+        t_interior = 2:hours_per_subperiod
+
+        if !isempty(stor_long_duration)
+            SOC_t[stor_long_duration, t_r] = socw[stor_long_duration, r] .* (1 .- self_discharge.(gen[stor_long_duration])) .+ efficiency_up.(gen[stor_long_duration]) .* v_charge[stor_long_duration, t_start_w] .- 1 ./ efficiency_down.(gen[stor_long_duration]) .* v_P[stor_long_duration, t_start_w]
+
+            for t_int in t_interior
+                t = hours_per_subperiod * (w - 1) + t_int
+                SOC_t[stor_long_duration, t_r + t_int - 1] = SOC_t[stor_long_duration, t_r + t_int - 2] .* (1 .- self_discharge.(gen[stor_long_duration])) .+ efficiency_up.(gen[stor_long_duration]) .* v_charge[stor_long_duration, t] .- 1 ./ efficiency_down.(gen[stor_long_duration]) .* v_P[stor_long_duration, t]
+            end
+        end
+
+        if !isempty(stor_hydro_long_duration)
+            SOC_t[stor_hydro_long_duration, t_r] = socw[stor_hydro_long_duration, r] .- 1 ./ efficiency_down.(gen[stor_long_duration]) .* v_P[stor_hydro_long_duration, t_start_w] .- v_spill[stor_hydro_long_duration, t_start_w] .+ pP_max[stor_hydro_long_duration, t_start_w] .* e_total_cap[stor_hydro_long_duration]
+
+            for t_int in t_interior
+                t = hours_per_subperiod * (w - 1) + t_int
+                SOC_t[stor_hydro_long_duration, t_r + t_int - 1] = SOC_t[stor_hydro_long_duration, t_r + t_int - 2] .- 1 ./ efficiency_down.(gen[stor_long_duration]) .* v_P[stor_hydro_long_duration, t] .- v_spill[stor_hydro_long_duration, t] .+ pP_max[stor_hydro_long_duration, t] .* e_total_cap[stor_hydro_long_duration]
+            end
+        end
+    end
+    df_SOC_t = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
+    df_SOC_t = hcat(df_SOC_t, DataFrame(SOC_t, :auto))
+    auxNew_Names = [Symbol("Resource"); Symbol("Zone"); [Symbol("n$t") for t in 1:T_hor]]
+    rename!(df_SOC_t,auxNew_Names)
+    CSV.write(joinpath(path, "StorageEvol.csv"), dftranspose(df_SOC_t, false), writeheader=false)
+
 end