From 3220650e396fffbe00f170211981c4db8bdaaa34 Mon Sep 17 00:00:00 2001 From: Alinson S Xavier Date: Thu, 20 Jan 2022 10:18:19 -0600 Subject: [PATCH] Implement new reserves --- docs/format.md | 2 +- docs/model.md | 2 +- instances/test/case14.json.gz | Bin 1834 -> 1842 bytes src/instance/read.jl | 11 +- src/instance/structs.jl | 2 + src/model/formulations/ArrCon2000/ramp.jl | 8 +- .../formulations/DamKucRajAta2016/ramp.jl | 10 +- src/model/formulations/Gar1962/prod.jl | 4 +- src/model/formulations/MorLatRam2013/ramp.jl | 12 +-- src/model/formulations/PanGua2016/ramp.jl | 8 +- src/model/formulations/base/system.jl | 24 +++++ src/model/formulations/base/unit.jl | 40 +++++-- src/solution/solution.jl | 14 +++ src/validation/validate.jl | 33 +++++- test/model/formulations_test.jl | 100 +++++++++++------- test/runtests.jl | 10 +- test/usage.jl | 2 +- 17 files changed, 201 insertions(+), 81 deletions(-) diff --git a/docs/format.md b/docs/format.md index c05b184..02ef89f 100644 --- a/docs/format.md +++ b/docs/format.md @@ -28,7 +28,7 @@ Each section is described in detail below. For a complete example, see [case14]( ### Parameters -This section describes system-wide parameters, such as power balance and reserve shortfall penalties, and optimization parameters, such as the length of the planning horizon and the time. +This section describes system-wide parameters, such as power balance penalty, and optimization parameters, such as the length of the planning horizon and the time. | Key | Description | Default | Time series? | :----------------------------- | :------------------------------------------------ | :------: | :------------: diff --git a/docs/model.md b/docs/model.md index 218dcd5..8936207 100644 --- a/docs/model.md +++ b/docs/model.md @@ -23,7 +23,7 @@ Name | Symbol | Description | Unit `switch_off[g,t]` | $w_{g}(t)$ | True if generator `g` switches off at time `t`. | Binary `prod_above[g,t]` |$p'_{g}(t)$ | Amount of power produced by generator `g` above its minimum power output at time `t`. For example, if the minimum power of generator `g` is 100 MW and `g` is producing 115 MW of power at time `t`, then `prod_above[g,t]` equals `15.0`. | MW `segprod[g,t,k]` | $p^k_g(t)$ | Amount of power from piecewise linear segment `k` produced by generator `g` at time `t`. For example, if cost curve for generator `g` is defined by the points `(100, 1400)`, `(110, 1600)`, `(130, 2200)` and `(135, 2400)`, and if the generator is producing 115 MW of power at time `t`, then `segprod[g,t,:]` equals `[10.0, 5.0, 0.0]`.| MW -`reserve[g,t]` | $r_g(t)$ | Amount of reserves provided by generator `g` at time `t`. | MW +`reserve[r,g,t]` | $r_g(t)$ | Amount of reserve `r` provided by unit `g` at time `t`. | MW `startup[g,t,s]` | $\delta^s_g(t)$ | True if generator `g` switches on at time `t` incurring start-up costs from start-up category `s`. | Binary diff --git a/instances/test/case14.json.gz b/instances/test/case14.json.gz index e4d49b8dc58b35c39d0d90ee180b9d33512ecdb4..682622f6c4d6585d5eff46687fdca625fd942c31 100644 GIT binary patch literal 1842 zcmV-22hI2&iwFoQi0NSf17l%xWid1^YIARH0PR{!Z`(!?zUNm6LJle5Y7|Le8y z+D>SGq2?#s8*i)fvRyP8zGG3=&1}(D@6*|<)8D)I{jnw6eUy|Z!$n0c0=a0n*8b$7 zc{u{ckZ2sWsx`Ta8g>FYYm%yISF_7BPj2h`Yaoyis)tL!dR2Gudi=~%N%mQL^i&1@~pWfuYF#r`S#_~kq+aL4xc~L0fuYN5s^N@ z>(7pKr4$Hn2?_ugg(M;!LR7-?BJzMpA1;L>rW-x74hNEgA7~pR1wxV}1xk)21w!z& zMVF}B;#=Hqsy5|thT?I?Ky)Gvre{ns+N++Y%PcXdH@S5KcCqSV-x zD60t{)%K($JM3J{1LiFW~#WmkfJcSDV$BRHx|6xMkPP7umQ zL)4nsMQJda1CwBwiR%YpwAa+_YvHiHCh*{n)P{5mC4z<+I8gn8!59vFF79g@umYT; zZVE6|t%*daf)dUv@)-;m0~R=@QF8BVl)@uJP%{{b@3odd6oqrz5P3cdE0s|!q!W}n zL~4*yoFu#oL`0m_y0_j@<2~V`iID|!Zw5!AQAjmJVI+cghTK4xpq4RGf@PuHpS+Qv zgCS$o3XQBs5VXP!6u}V=fjin_Zen5Go4b`40oj(gNIZp)&Qe&dh!NS6-0>jCv2N7@ z=T%bF%dD=u@+Z%VeqoRkUVK)z)gql;j1BV7lj=HcHv5Aq>$@~rG)aLfUX?e?`ktC> zR9-%Q`rYYSTQAbp{bOg{GaA!rbiqmQ0pmmV^Yr|;vRYzsb(NPNw)<=tP7LNhF0tQD zi^VM=ZwUFH9;H@kk>m~b$}hjZlULiLho|0=6XCX}JuSjP-(Y0yNyw~FR z)p9?(I}x=7CuPL8PQwTrQl)9{00V%%1E@ik8X@rz8P4Lp!8TsoE?G{iO?ohxe>WkV@x=mHYZW@vOR+6j}FL?#6DSKv+s21a#k40J8C*} z4y}181M2?3o^);*uw@79wZoq=Cz;NdLxww$N@8Txgvi{WLcnNL8vCyiFVWNaa)|0M zHvtnEosR8B&=SikjR6UupZqD2uci zAqyhcv9xQJb|NFOv;#{!kvCY{4NE)0?=0<>rJbN*mUhR|PS7PwyJu-9NR6dEu(T80 z!_ppE+NqkF7fzC=|KM=t#Vb{u^7Qi}mMS87`gt)+l}$YTyttz38=ih%WKop}Pd_iV zsEUE7pBG(JYUkc?o4azP3{3+bUb6uj{m^vnG3=&hj$3?4`_Az2ncpiL&H< zTD0lbRCO?g{UKc@#b&Z@x2CRJ@b14vueZar`^tQP==^q-c970iSy5p2?*n_YEZbr; g8>Ud5-IP^xmE`&KX!s7~i%&290X+DJi-Ru!0L~KgVw9v6(CE79!a-B{(DFJwm7oJPB!hKPA*GQqmRS+zL}x@;b?)M<(sU@ z)_IpV?ed$&hwHKH@snbmFD~n*_`9wa&o95iS#v%6@x{%9Kj&56WL-Ue@Ir0)(g>{js%n_$Vn~q=SkY1ai=Bjs3|( z^K1f)A<;N$RBLh+E$jqz++|JIZx&~HnO(KF$3P$=gnP+JedCY;!z0d=@8VP7=pwmjZep-A)I7( zdt$%vBcseZK?R!e00Ib@9-i5;`$~(xj>n{6_W$<-T73+SDD7!&#Z*_KJ(O_}cvh1&$Wm!>PEe3VIxEfx4y}0Vj za_9PX*jh`kq$wI$1ReW2Z}TCl^Rl=oPK&bWuE=9wHfp(jxHQvYn(6TQnGTSSIY%&k zM%O=?=}IXO-Vzc3E(%F79YR#X@)&sl)6=2wgy}|4jKhJV;1g?8qCl)9MS)sJiUP6l zb&ny@^wqaG-8FsA(G11YjAYpt)(rPw&%;%9zsK+t zu6Ix8`p@iAd9C7KrfWEQ+I~_l%U;xO0*!4pG&u`v2<5f6A`10 zFitC?rPCmE^qLAkEusw3GxYG3sG)X4RSKh!g6*n<*ya7JoFhJ_MALkt|K{y;K@1D}igx(2YDc8)qJY{0cO zk)SFl;k+WBLBbfYz;TU|dmp0|9uB6w%W2@DBp7?~0*3+4XejRGAE$x$mb zsvaTG3NugyM+5}UXp6atg>`T3R$c^DTjC&b7d|>mVYMPgR7-NkgPg{?RSTRnS=FwK zwjHXLvZ%&|K}kgMaosnod~rH;+kKKX7kPI-pIk-0&9haPRcPXQeYtLLrO8I+<-=#( z!X5YRD&O2bc04?zF`Y&iob(=$A99?hm%r7`8f&8SvVOnaj6&L0FMnKPQ|x7S^m@-qlLT!;V(E=cGT%b2V zFmbeje?f$vr_E9%%1i-qTTp7*HuAMCGE6 zng&f2R(pcxg^I$OYkbR@g2rqq!CNZ)cww?w})O&k^_B<_-QK|Uw;m{H#L^|v6 zzlug0dILn^x#tEVYXcTh%ttt2-xd(GaX}8D?k3z$)u1ETNid9HCl_$E->sJW`Q3@A zEjTG7wsjgN)Q~DoM;k~0jy9l@EHxqGAu@u+dxLGfwgW%0Hk&{_pE_qz^J)r_EW^yljsk##aXu5hZ#)U#6%I3llJb(dpRX z3R+@WrIC;T`qA&Dw!7o+>#8fNp)D-Njq!@ep|8%*^HrCh-8a(2mv<%W%exhk>sZz` z%Q}&fSk{4MoyZ$3>xN~W(07(~%d$?$Fw44QStsO@W!J+0&pwY^RK>uv&tn&r z+IjYQ{GxI&FQH8HYb#a0X^K_;vdycu=!$pwqO7yCQOex3JAN}wlqK)-s?WD>Rg+uT zAM$lp-QTR+ty|YEc(*^UX}4)Q{FEj^baJ)H2S~@8qN*_akAb~f*L`(=H%y^AzO0+> YJS)rjtKl1vM<0*=0UC%}#;7g;0NAF7&j0`b diff --git a/src/instance/read.jl b/src/instance/read.jl index a7c9359..5606adb 100644 --- a/src/instance/read.jl +++ b/src/instance/read.jl @@ -125,6 +125,11 @@ function _from_json(json; repair = true) name = reserve_name, type = lowercase(dict["Type"]), amount = timeseries(dict["Amount (MW)"]), + units = [], + shortfall_penalty = scalar( + dict["Shortfall penalty (\$/MW)"], + default = -1, + ), ) name_to_reserve[reserve_name] = reserve push!(reserves2, reserve) @@ -171,7 +176,8 @@ function _from_json(json; repair = true) # Read reserves unit_reserves = Reserve[] if "Reserve eligibility" in keys(dict) - unit_reserves = [name_to_reserve[n] for n in dict["Reserve eligibility"]] + unit_reserves = + [name_to_reserve[n] for n in dict["Reserve eligibility"]] end # Read and validate initial conditions @@ -215,6 +221,9 @@ function _from_json(json; repair = true) unit_reserves, ) push!(bus.units, unit) + for r in unit_reserves + push!(r.units, unit) + end name_to_unit[unit_name] = unit push!(units, unit) end diff --git a/src/instance/structs.jl b/src/instance/structs.jl index 8b92777..33827ba 100644 --- a/src/instance/structs.jl +++ b/src/instance/structs.jl @@ -24,6 +24,8 @@ Base.@kwdef mutable struct Reserve name::String type::String amount::Vector{Float64} + units::Vector + shortfall_penalty::Float64 end mutable struct Unit diff --git a/src/model/formulations/ArrCon2000/ramp.jl b/src/model/formulations/ArrCon2000/ramp.jl index e9bb288..080abe6 100644 --- a/src/model/formulations/ArrCon2000/ramp.jl +++ b/src/model/formulations/ArrCon2000/ramp.jl @@ -19,10 +19,10 @@ function _add_ramp_eqs!( RD = g.ramp_down_limit SU = g.startup_limit SD = g.shutdown_limit - reserve = model[:reserve] eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_ramp_up) is_initially_on = (g.initial_status > 0) + reserve = _total_reserves(model, g) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -41,7 +41,7 @@ function _add_ramp_eqs!( model, g.min_power[t] + prod_above[gn, t] + - (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <= + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) end @@ -51,7 +51,7 @@ function _add_ramp_eqs!( prod_above[gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? - reserve[gn, t] : 0.0 + reserve[t] : 0.0 ) min_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] @@ -82,7 +82,7 @@ function _add_ramp_eqs!( prod_above[gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? - reserve[gn, t-1] : 0.0 + reserve[t-1] : 0.0 ) min_prod_this_period = g.min_power[t] * is_on[gn, t] + prod_above[gn, t] diff --git a/src/model/formulations/DamKucRajAta2016/ramp.jl b/src/model/formulations/DamKucRajAta2016/ramp.jl index 9afd247..cee2db9 100644 --- a/src/model/formulations/DamKucRajAta2016/ramp.jl +++ b/src/model/formulations/DamKucRajAta2016/ramp.jl @@ -23,7 +23,7 @@ function _add_ramp_eqs!( gn = g.name eq_str_ramp_down = _init(model, :eq_str_ramp_down) eq_str_ramp_up = _init(model, :eq_str_ramp_up) - reserve = model[:reserve] + reserve = _total_reserves(model, g) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -48,10 +48,8 @@ function _add_ramp_eqs!( # end max_prod_this_period = - prod_above[gn, t] + ( - RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? - reserve[gn, t] : 0.0 - ) + prod_above[gn, t] + + (RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) min_prod_last_period = 0.0 if t > 1 && time_invariant min_prod_last_period = prod_above[gn, t-1] @@ -88,7 +86,7 @@ function _add_ramp_eqs!( max_prod_last_period = min_prod_last_period + ( t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ? - reserve[gn, t-1] : 0.0 + reserve[t-1] : 0.0 ) min_prod_this_period = prod_above[gn, t] on_last_period = 0.0 diff --git a/src/model/formulations/Gar1962/prod.jl b/src/model/formulations/Gar1962/prod.jl index e39a90e..fa83e79 100644 --- a/src/model/formulations/Gar1962/prod.jl +++ b/src/model/formulations/Gar1962/prod.jl @@ -26,7 +26,7 @@ function _add_production_limit_eqs!( eq_prod_limit = _init(model, :eq_prod_limit) is_on = model[:is_on] prod_above = model[:prod_above] - reserve = model[:reserve] + reserve = _total_reserves(model, g) gn = g.name for t in 1:model[:instance].time # Objective function terms for production costs @@ -44,7 +44,7 @@ function _add_production_limit_eqs!( end eq_prod_limit[gn, t] = @constraint( model, - prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t] + prod_above[gn, t] + reserve[t] <= power_diff * is_on[gn, t] ) end end diff --git a/src/model/formulations/MorLatRam2013/ramp.jl b/src/model/formulations/MorLatRam2013/ramp.jl index cbb8f94..5eaa178 100644 --- a/src/model/formulations/MorLatRam2013/ramp.jl +++ b/src/model/formulations/MorLatRam2013/ramp.jl @@ -22,7 +22,7 @@ function _add_ramp_eqs!( gn = g.name eq_ramp_down = _init(model, :eq_ramp_down) eq_ramp_up = _init(model, :eq_str_ramp_up) - reserve = model[:reserve] + reserve = _total_reserves(model, g) # Gar1962.ProdVars prod_above = model[:prod_above] @@ -43,7 +43,7 @@ function _add_ramp_eqs!( model, g.min_power[t] + prod_above[gn, t] + - (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <= + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <= g.initial_power + RU ) end @@ -61,7 +61,7 @@ function _add_ramp_eqs!( prod_above[gn, t] + ( RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? - reserve[gn, t] : 0.0 + reserve[t] : 0.0 ) min_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1] @@ -77,7 +77,7 @@ function _add_ramp_eqs!( eq_ramp_up[gn, t] = @constraint( model, prod_above[gn, t] + - (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) - + (RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) - prod_above[gn, t-1] <= RU ) end @@ -105,7 +105,7 @@ function _add_ramp_eqs!( prod_above[gn, t-1] + ( RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? - reserve[gn, t-1] : 0.0 + reserve[t-1] : 0.0 ) min_prod_this_period = g.min_power[t] * is_on[gn, t] + prod_above[gn, t] @@ -121,7 +121,7 @@ function _add_ramp_eqs!( eq_ramp_down[gn, t] = @constraint( model, prod_above[gn, t-1] + - (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t-1] : 0.0) - + (RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) - prod_above[gn, t] <= RD ) end diff --git a/src/model/formulations/PanGua2016/ramp.jl b/src/model/formulations/PanGua2016/ramp.jl index f040824..fe5d334 100644 --- a/src/model/formulations/PanGua2016/ramp.jl +++ b/src/model/formulations/PanGua2016/ramp.jl @@ -12,7 +12,7 @@ function _add_ramp_eqs!( # TODO: Move upper case constants to model[:instance] RESERVES_WHEN_SHUT_DOWN = true gn = g.name - reserve = model[:reserve] + reserve = _total_reserves(model, g) eq_str_prod_limit = _init(model, :eq_str_prod_limit) eq_prod_limit_ramp_up_extra_period = _init(model, :eq_prod_limit_ramp_up_extra_period) @@ -56,7 +56,7 @@ function _add_ramp_eqs!( model, prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + - reserve[gn, t] <= + reserve[t] <= Pbar * is_on[gn, t] - (t < T ? (Pbar - SD) * switch_off[gn, t+1] : 0.0) - sum( (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for @@ -71,7 +71,7 @@ function _add_ramp_eqs!( model, prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + - reserve[gn, t] <= + reserve[t] <= Pbar * is_on[gn, t] - sum( (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for i in 0:min(UT - 1, TRU, t - 1) @@ -88,7 +88,7 @@ function _add_ramp_eqs!( model, prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + - (RESERVES_WHEN_SHUT_DOWN ? reserve[gn, t] : 0.0) <= + (RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <= Pbar * is_on[gn, t] - sum( (Pbar - (SD + i * RD)) * switch_off[gn, t+1+i] for i in 0:KSD diff --git a/src/model/formulations/base/system.jl b/src/model/formulations/base/system.jl index d6bf573..7fef70f 100644 --- a/src/model/formulations/base/system.jl +++ b/src/model/formulations/base/system.jl @@ -52,5 +52,29 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing ) end end + + eq_min_reserve2 = _init(model, :eq_min_reserve2) + for r in instance.reserves2 + for t in 1:instance.time + # Equation (68) in Kneuven et al. (2020) + # As in Morales-España et al. (2013a) + # Akin to the alternative formulation with max_power_avail + # from Carrión and Arroyo (2006) and Ostrowski et al. (2012) + eq_min_reserve2[r.name, t] = @constraint( + model, + sum(model[:reserve2][r.name, g.name, t] for g in r.units) + + model[:reserve_shortfall2][r.name, t] >= r.amount[t] + ) + + # Account for shortfall contribution to objective + if r.shortfall_penalty >= 0 + add_to_expression!( + model[:obj], + r.shortfall_penalty, + model[:reserve_shortfall2][r.name, t], + ) + end + end + end return end diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl index e701977..51de4f2 100644 --- a/src/model/formulations/base/unit.jl +++ b/src/model/formulations/base/unit.jl @@ -55,13 +55,19 @@ function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing (model[:instance].shortfall_penalty[t] >= 0) ? @variable(model, lower_bound = 0) : 0.0 end - return -end -function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing - reserve = model[:reserve] - for t in 1:model[:instance].time - add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0) + reserve2 = _init(model, :reserve2) + reserve_shortfall2 = _init(model, :reserve_shortfall2) + for r in g.reserves + for t in 1:model[:instance].time + reserve2[r.name, g.name, t] = @variable(model, lower_bound = 0) + if (r.name, t) ∉ keys(reserve_shortfall2) + reserve_shortfall2[r.name, t] = @variable(model, lower_bound = 0) + if r.shortfall_penalty < 0 + set_upper_bound(reserve_shortfall2[r.name, t], 0.0) + end + end + end end return end @@ -81,7 +87,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing eq_startup_limit = _init(model, :eq_startup_limit) is_on = model[:is_on] prod_above = model[:prod_above] - reserve = model[:reserve] + reserve = _total_reserves(model, g) switch_off = model[:switch_off] switch_on = model[:switch_on] T = model[:instance].time @@ -89,7 +95,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing # Startup limit eq_startup_limit[g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[g.name, t] <= + prod_above[g.name, t] + reserve[t] <= (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] - max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t] ) @@ -117,7 +123,7 @@ function _add_ramp_eqs!( formulation::RampingFormulation, )::Nothing prod_above = model[:prod_above] - reserve = model[:reserve] + reserve = _total_reserves(model, g) eq_ramp_up = _init(model, :eq_ramp_up) eq_ramp_down = _init(model, :eq_ramp_down) for t in 1:model[:instance].time @@ -126,14 +132,14 @@ function _add_ramp_eqs!( if _is_initially_on(g) == 1 eq_ramp_up[g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[g.name, t] <= + prod_above[g.name, t] + reserve[t] <= (g.initial_power - g.min_power[t]) + g.ramp_up_limit ) end else eq_ramp_up[g.name, t] = @constraint( model, - prod_above[g.name, t] + reserve[g.name, t] <= + prod_above[g.name, t] + reserve[t] <= prod_above[g.name, t-1] + g.ramp_up_limit ) end @@ -216,3 +222,15 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing ) end end + +function _total_reserves(model, g)::Vector + T = model[:instance].time + reserve = [model[:reserve][g.name, t] for t in 1:T] + if !isempty(g.reserves) + reserve += [ + sum(model[:reserve2][r.name, g.name, t] for r in g.reserves) for + t in 1:model[:instance].time + ] + end + return reserve +end diff --git a/src/solution/solution.jl b/src/solution/solution.jl index 5fd8bbd..35de2b3 100644 --- a/src/solution/solution.jl +++ b/src/solution/solution.jl @@ -67,5 +67,19 @@ function solution(model::JuMP.Model)::OrderedDict sol["Price-sensitive loads (MW)"] = timeseries(model[:loads], instance.price_sensitive_loads) end + sol["Reserve 2 (MW)"] = OrderedDict( + r.name => OrderedDict( + g.name => [ + value(model[:reserve2][r.name, g.name, t]) for + t in 1:instance.time + ] for g in r.units + ) for r in instance.reserves2 + ) + sol["Reserve shortfall 2 (MW)"] = OrderedDict( + r.name => [ + value(model[:reserve_shortfall2][r.name, t]) for + t in 1:instance.time + ] for r in instance.reserves2 + ) return sol end diff --git a/src/validation/validate.jl b/src/validation/validate.jl index d342ef9..62def52 100644 --- a/src/validation/validate.jl +++ b/src/validation/validate.jl @@ -46,6 +46,12 @@ function _validate_units(instance, solution; tol = 0.01) for unit in instance.units production = solution["Production (MW)"][unit.name] reserve = solution["Reserve (MW)"][unit.name] + if !isempty(unit.reserves) + reserve += sum( + solution["Reserve 2 (MW)"][r.name][unit.name] for + r in unit.reserves + ) + end actual_production_cost = solution["Production cost (\$)"][unit.name] actual_startup_cost = solution["Startup cost (\$)"][unit.name] is_on = bin(solution["Is on"][unit.name]) @@ -137,9 +143,11 @@ function _validate_units(instance, solution; tol = 0.01) # If unit is off, must produce zero if !is_on[t] && production[t] + reserve[t] > tol @error @sprintf( - "Unit %s produces power at time %d while off", + "Unit %s produces power at time %d while off (%.2f + %.2f > 0)", unit.name, - t + t, + production[t], + reserve[t], ) err_count += 1 end @@ -338,6 +346,27 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01) ) err_count += 1 end + + # Verify reserves + for r in instance.reserves2 + provided = sum( + solution["Reserve 2 (MW)"][r.name][g.name][t] for g in r.units + ) + shortfall = solution["Reserve shortfall 2 (MW)"][r.name][t] + required = r.amount[t] + + if provided + shortfall < required - tol + @error @sprintf( + "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)", + r.name, + t, + provided, + shortfall, + required, + ) + err_count += 1 + end + end end return err_count diff --git a/test/model/formulations_test.jl b/test/model/formulations_test.jl index 3b08dc5..eba99e0 100644 --- a/test/model/formulations_test.jl +++ b/test/model/formulations_test.jl @@ -4,6 +4,8 @@ using UnitCommitment using JuMP +using Cbc +using JSON import UnitCommitment: ArrCon2000, CarArr2006, @@ -19,56 +21,78 @@ if ENABLE_LARGE_TESTS using Gurobi end -function _small_test(formulation::Formulation)::Nothing - instances = ["matpower/case118/2017-02-01", "test/case14"] - for instance in instances - # Should not crash - UnitCommitment.build_model( - instance = UnitCommitment.read_benchmark(instance), - formulation = formulation, - ) +function _small_test(formulation::Formulation; dump::Bool = false)::Nothing + instance = UnitCommitment.read_benchmark("test/case14") + model = UnitCommitment.build_model( + instance = instance, + formulation = formulation, + optimizer = Cbc.Optimizer, + variable_names = true, + ) + UnitCommitment.optimize!(model) + solution = UnitCommitment.solution(model) + if dump + open("/tmp/ucjl.json", "w") do f + return write(f, JSON.json(solution, 2)) + end + write_to_file(model, "/tmp/ucjl.lp") end + @test UnitCommitment.validate(instance, solution) return end function _large_test(formulation::Formulation)::Nothing - instances = ["pglib-uc/ca/Scenario400_reserves_1"] - for instance in instances - instance = UnitCommitment.read_benchmark(instance) - model = UnitCommitment.build_model( - instance = instance, - formulation = formulation, - optimizer = Gurobi.Optimizer, - ) - UnitCommitment.optimize!( - model, - XavQiuWanThi2019.Method(two_phase_gap = false, gap_limit = 0.1), - ) - solution = UnitCommitment.solution(model) - @test UnitCommitment.validate(instance, solution) - end + instance = + UnitCommitment.read_benchmark("pglib-uc/ca/Scenario400_reserves_1") + model = UnitCommitment.build_model( + instance = instance, + formulation = formulation, + optimizer = Gurobi.Optimizer, + ) + UnitCommitment.optimize!( + model, + XavQiuWanThi2019.Method(two_phase_gap = false, gap_limit = 0.1), + ) + solution = UnitCommitment.solution(model) + @test UnitCommitment.validate(instance, solution) return end -function _test(formulation::Formulation)::Nothing - _small_test(formulation) +function _test(formulation::Formulation; dump::Bool = false)::Nothing + _small_test(formulation; dump) if ENABLE_LARGE_TESTS _large_test(formulation) end end @testset "formulations" begin - _test(Formulation()) - _test(Formulation(ramping = ArrCon2000.Ramping())) - # _test(Formulation(ramping = DamKucRajAta2016.Ramping())) - _test( - Formulation( - ramping = MorLatRam2013.Ramping(), - startup_costs = MorLatRam2013.StartupCosts(), - ), - ) - _test(Formulation(ramping = PanGua2016.Ramping())) - _test(Formulation(pwl_costs = Gar1962.PwlCosts())) - _test(Formulation(pwl_costs = CarArr2006.PwlCosts())) - _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts())) + @testset "default" begin + _test(Formulation()) + end + @testset "ArrCon2000" begin + _test(Formulation(ramping = ArrCon2000.Ramping())) + end + @testset "DamKucRajAta2016" begin + _test(Formulation(ramping = DamKucRajAta2016.Ramping())) + end + @testset "MorLatRam2013" begin + _test( + Formulation( + ramping = MorLatRam2013.Ramping(), + startup_costs = MorLatRam2013.StartupCosts(), + ), + ) + end + @testset "PanGua2016" begin + _test(Formulation(ramping = PanGua2016.Ramping())) + end + @testset "Gar1962" begin + _test(Formulation(pwl_costs = Gar1962.PwlCosts())) + end + @testset "CarArr2006" begin + _test(Formulation(pwl_costs = CarArr2006.PwlCosts())) + end + @testset "KnuOstWat2018" begin + _test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts())) + end end diff --git a/test/runtests.jl b/test/runtests.jl index 836eda2..7f1a1be 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,10 +21,12 @@ const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV)) @testset "model" begin include("model/formulations_test.jl") end - @testset "XavQiuWanThi19" begin - include("solution/methods/XavQiuWanThi19/filter_test.jl") - include("solution/methods/XavQiuWanThi19/find_test.jl") - include("solution/methods/XavQiuWanThi19/sensitivity_test.jl") + @testset "solution" begin + @testset "XavQiuWanThi19" begin + include("solution/methods/XavQiuWanThi19/filter_test.jl") + include("solution/methods/XavQiuWanThi19/find_test.jl") + include("solution/methods/XavQiuWanThi19/sensitivity_test.jl") + end end @testset "transform" begin include("transform/initcond_test.jl") diff --git a/test/usage.jl b/test/usage.jl index 5ff589d..0cbf799 100644 --- a/test/usage.jl +++ b/test/usage.jl @@ -4,7 +4,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON -@testset "build_model" begin +@testset "usage" begin instance = UnitCommitment.read_benchmark("test/case14") for line in instance.lines, t in 1:4 line.normal_flow_limit[t] = 10.0