Merge pull request #36 from NVE/elasticdemand

Elasticdemand

src/TuLiPa.jl

@@ -43,10 +43,16 @@ include("trait_metadata.jl")
 include("parameters.jl")
 
 # Toplevel model objects
+
+
+
 include("obj_balance.jl")
 include("obj_flow.jl")
 include("obj_storage.jl")
 include("obj_aggsupplycurve.jl")
+
+include("obj_elastic_demand.jl")
+
 include("trait_softbound.jl")
 include("trait_startupcost.jl")
 include("trait_ramping.jl")

src/abstracttypes.jl

@@ -546,7 +546,18 @@ abstract type TimeVector end
 
 abstract type Offset end
 
-
-
-
-
+# ---- ElasticDemand -------------------
+# 
+# Is used to create demand that should change with price
+# 
+# Interface:
+#    build!(prob, demand)
+#    setconstants!(prob, demand)
+#    update!(prob, demand, start)
+#    assemble!(demand)
+#    getid(demand)
+#    getbalance(demand)
+#    getdemand(prob, obj, timeix)
+#    getparent(var::BaseElasticDemand)
+
+abstract type ElasticDemand end

src/input_utils.jl

@@ -61,6 +61,7 @@ const CAPACITY_CONCEPT = "Capacity"
 const PRICE_CONCEPT = "Price"
 const CONVERSION_CONCEPT = "Conversion"
 const LOSS_CONCEPT = "Loss"
+const ELASTIC_DEMAND_CONCEPT = "ElasticDemand"
 
 # ---- Functions to parse input files containing data elements ----
 

src/obj_elastic_demand.jl

@@ -0,0 +1,243 @@
+mutable struct BaseElasticDemand{B, P} <: ElasticDemand
+    id::Id
+    balance::B
+    firm_demand::P
+    normal_price::Float64
+    price_elasticity::Float64
+    min_price::Float64 
+    max_price::Float64
+    N::Int64
+    segment_capacities::Vector{Float64}
+    reserve_prices::Vector{Float64}
+
+    function BaseElasticDemand(
+            id::Id, 
+            balance::B, 
+            firm_demand_param::P, 
+            price_elasticity::Float64, 
+            normal_price::Float64, 
+            max_price::Float64, 
+            min_price::Float64,
+            threshold::Float64
+        ) where {B <: Balance, P <: Param} 
+
+        min_relative_demand = price_to_relative_demand(normal_price, price_elasticity, max_price)
+        max_relative_demand = price_to_relative_demand(normal_price, price_elasticity, min_price)
+
+        L, reserve_prices, N = optimize_segments(
+            normal_price, 
+            price_elasticity, 
+            min_relative_demand, 
+            max_relative_demand, 
+            threshold
+        ) 
+
+        reserve_prices = adjust_prices_for_demand_curve_area(
+            normal_price, 
+            price_elasticity, 
+            min_relative_demand, 
+            max_relative_demand,
+            L, 
+            reserve_prices
+        )
+
+        segment_capacities = [first(L), diff(L)...]
+
+        new{B, P}(
+            id, 
+            balance, 
+            firm_demand_param, 
+            normal_price, 
+            price_elasticity, 
+            min_price, 
+            max_price, 
+            N,
+            segment_capacities, 
+            reserve_prices
+        )
+    end
+end
+
+getid(var::BaseElasticDemand) = var.id
+getbalance(var::BaseElasticDemand) = var.balance
+getparent(var::BaseElasticDemand) = var.balance
+
+function getdemand(p::Prob, var::BaseElasticDemand, timeix::Int64)
+    total_demand = 0
+    for seg_no in 1:var.N
+        seg_id = create_segment_id(var, seg_no)
+        total_demand += getvarvalue(prob, seg_id, timeix)
+    end
+    return total_demand
+end
+
+# Expects f divided by f_ref
+function relative_demand_to_price(normal_price, price_elasticity, f)
+    p_ref = normal_price
+    e = price_elasticity
+    return p_ref .* f.^(1/e)
+end
+
+# Price_to_demand divided by f_ref
+function price_to_relative_demand(normal_price, price_elasticity, p)
+    p_ref = normal_price
+    e = price_elasticity
+    return (p./p_ref).^e
+end
+
+function create_segment_id(var::BaseElasticDemand, seg_no::Int)
+    return Id(var.id.conceptname, string(var.id.instancename, seg_no))
+end
+
+function adjust_prices_for_demand_curve_area(
+        normal_price, 
+        price_elasticity, 
+        min_relative_demand, 
+        max_relative_demand,
+        L, 
+        reserve_prices
+    )
+    max_price = reserve_prices[1]
+    demand_integral = (p_normal, e, f) -> (p_normal .* e .* f.^(1 + 1/e))/(1 + e)
+    demand_area = (p_normal, e, f1, f2) -> demand_integral(p_normal, e, f2) - demand_integral(p_normal, e, f1)
+    demand_approx_area = (p, seg) -> sum(p .* seg)
+    factor = (p_normal, e, f1, f2, p, seg) -> demand_area(p_normal, e, f1, f2) / demand_approx_area(p, seg)
+    reserve = reserve_prices[1:end-1]
+    seg = diff(L)
+    p_factor = factor(
+        normal_price, 
+        price_elasticity, 
+        min_relative_demand, 
+        max_relative_demand, 
+        reserve, 
+        seg
+    ) 
+    reserve_prices = reserve * p_factor
+    insert!(reserve_prices, 1, max_price)
+    return reserve_prices
+end
+
+function optimize_segments(normal_price, price_elasticity, min_relative_demand, max_relative_demand, tolerance; max_depth = 6)
+    f = (x) -> relative_demand_to_price(normal_price, price_elasticity, x)
+    x_points = adaptive_sampling(f, min_relative_demand, max_relative_demand, tolerance, max_depth)
+    y_points = relative_demand_to_price(normal_price, price_elasticity, x_points)
+    N = length(x_points)
+    @assert N <= 10
+    return x_points, y_points, N
+end
+
+function adaptive_sampling(f, x_start, x_end, tolerance, max_depth)
+    x_points = [x_start, x_end]
+    y_points = [f(x_start), f(x_end)]
+    tolerance *= abs(y_points[1] - y_points[2])
+    function recursive_sampling(x0, x1, y0, y1, depth)
+        depth > max_depth && return
+        x_mid = (x0 + x1) / 2
+        y_mid = f(x_mid)
+        y_interp = (y0 + y1) / 2
+        if abs(y0 - y_interp) > tolerance
+            push!(x_points, x_mid)
+            push!(y_points, y_mid)
+            recursive_sampling(x0, x_mid, y0, y_mid, depth + 1)
+            recursive_sampling(x_mid, x1, y_mid, y1, depth + 1)
+        end
+    end
+    recursive_sampling(x_start, x_end, y_points[1], y_points[2], 0)
+    sort!(x_points)
+    return x_points
+end
+
+function build!(p::Prob, var::BaseElasticDemand)
+    T = getnumperiods(gethorizon(var.balance))
+    for i in 1:var.N
+        addvar!(p, create_segment_id(var, i), T)
+    end
+end
+
+function setconstants!(p::Prob, var::BaseElasticDemand)
+    balanceid = getid(getbalance(var))
+    T = getnumperiods(gethorizon(var.balance))
+    for i in 1:var.N
+        varid = create_segment_id(var, i)
+        for t in 1:T
+            setconcoeff!(p, balanceid, varid, t, t, 1.0)
+        end
+        varid = create_segment_id(var, i)
+        for t in 1:T
+            setobjcoeff!(p, varid, t, -var.reserve_prices[i])
+            setlb!(p, varid, t, 0.0)
+        end
+    end
+end
+
+function _update_ub(p, horizon, start, var, varid, t, i)
+    querystart = getstarttime(horizon, t, start)
+    querydelta = gettimedelta(horizon, t)
+    value = getparamvalue(var.firm_demand, querystart, querydelta)
+    setub!(p, varid, t, value * var.segment_capacities[i])
+end
+
+function stateful_update(p, horizon, start, var, varid, T, i)
+    for t in 1:T
+        _update_ub(p, horizon, start, var, varid, t, i)
+    end
+end
+
+function non_stateful_update(p, horizon, start, var, varid, T, i)
+    for t in 1:T
+        (future_t, ok) = mayshiftfrom(horizon, t)
+        if ok
+            value = getub(p, varid, future_t)
+            setub!(p, varid, t, value)
+        end
+    end
+    for t in 1:T
+        if mustupdate(horizon, t)
+            _update_ub(p, horizon, start, var, varid, t, i)
+        end
+    end
+end
+
+function update!(p::Prob, var::BaseElasticDemand, start::ProbTime)
+    horizon = gethorizon(var.balance)
+    T = getnumperiods(horizon)
+    update_func = isstateful(var.firm_demand) ? stateful_update : non_stateful_update
+    for i in 1:var.N
+        varid = create_segment_id(var, i)
+        update_func(p, horizon, start, var, varid, T, i)
+    end
+end
+
+function assemble!(var::BaseElasticDemand)::Bool   
+    return true
+end
+
+function includeBaseElasticDemand!(toplevel::Dict, lowlevel::Dict, elkey::ElementKey, value::Dict)
+    checkkey(toplevel, elkey)
+    deps = Id[]
+    balancename = getdictvalue(value, BALANCE_CONCEPT, String, elkey)
+    balancekey = Id(BALANCE_CONCEPT, balancename)
+    push!(deps, balancekey)
+
+    (id, firm_demand_param, ok) = getdictparamvalue(lowlevel, elkey, value)
+    push!(deps, id)
+
+    ok || return (false, deps)
+    haskey(toplevel, balancekey) || return (false, deps)
+
+    objkey = getobjkey(elkey)
+    price_elasticity = getdictvalue(value, "price_elasticity", Float64, elkey)
+    normal_price = getdictvalue(value, "normal_price", Float64, elkey)
+    max_price = getdictvalue(value, "max_price", Float64, elkey)
+    min_price = getdictvalue(value, "min_price", Float64, elkey)
+    threshold = getdictvalue(value, "threshold", Float64, elkey)
+
+    @assert min_price <= normal_price <= max_price
+
+    toplevel[objkey] = BaseElasticDemand(objkey, toplevel[balancekey], 
+        firm_demand_param, price_elasticity, normal_price, max_price, min_price, threshold
+    )
+    return (true, deps)    
+end
+
+INCLUDEELEMENT[TypeKey(ELASTIC_DEMAND_CONCEPT, "BaseElasticDemand")] = includeBaseElasticDemand!

src/reasoning_modelobjects.jl

@@ -480,6 +480,7 @@ replacebalance!(x::StartEqualStop, coupling, modelobjects) = nothing
 replacebalance!(x::BaseSoftBound, coupling, modelobjects) = nothing
 replacebalance!(x::HydroRampingWithout, coupling, modelobjects) = nothing
 replacebalance!(x::TransmissionRamping, coupling, modelobjects) = nothing # handled in replacebalance!(x::BaseFlow
+replacebalance!(x::BaseElasticDemand, coupling, modelobjects) = nothing
 
 function replacebalance!(x::BaseStorage, coupling, modelobjects)
     if haskey(coupling, getbalance(x))

src/results.jl

@@ -133,6 +133,15 @@ function order_result_objects(resultobjects, includeexogenprice=true)
                 push!(plantbalances,getid(balance))
             end
         end
+        if obj isa BaseElasticDemand
+            instance = getinstancename(getid(obj))
+            concept = getconceptname(getid(obj))
+            balance = getbalance(obj)
+            for c in 1:obj.N
+                push!(demands, create_segment_id(obj, c))
+                push!(demandbalances, getid(balance))
+            end
+        end
     end
     return powerbalances, rhsterms, rhstermbalances, plants, plantbalances, plantarrows, demands, demandbalances, demandarrows, hydrostorages, batterystorages
 end
@@ -219,16 +228,20 @@ function get_results!(problem, prices, rhstermvalues, production, consumption, h
 
         # Collect demand of all demands
         for i in 1:length(demands) # TODO: Balance and variable can have different horizons
-            if isexogen(modelobjects[demandbalances[i]])
-                arrow = demandarrows[demands[i]]
-                horizon = gethorizon(arrow)
-                conversionparam = getcontributionparam(arrow)
-                querytime = getstarttime(horizon, j, t)
-                querydelta = gettimedelta(horizon, j)
-                conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
-                consumption[jj, i] = getvarvalue(problem, demands[i], j)*conversionvalue/timefactor
+            if getconceptname(demands[i]) != ELASTIC_DEMAND_CONCEPT
+                if isexogen(modelobjects[demandbalances[i]])
+                    arrow = demandarrows[demands[i]]
+                    horizon = gethorizon(arrow)
+                    conversionparam = getcontributionparam(arrow)
+                    querytime = getstarttime(horizon, j, t)
+                    querydelta = gettimedelta(horizon, j)
+                    conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
+                    consumption[jj, i] = getvarvalue(problem, demands[i], j)*conversionvalue/timefactor
+                else
+                    consumption[jj, i] = getvarvalue(problem, demands[i], j)*abs(getconcoeff(problem, demandbalances[i], demands[i], j, j))/timefactor
+                end
             else
-                consumption[jj, i] = getvarvalue(problem, demands[i], j)*abs(getconcoeff(problem, demandbalances[i], demands[i], j, j))/timefactor
+                consumption[jj, i] = getvarvalue(problem, demands[i], j)/timefactor
             end
         end
 

test/includeelement_functions.jl

@@ -137,6 +137,8 @@ function test_all_includeelement_methods()
     test_includeBaseRHSTerm!()
     test_includeBaseSoftBound!()
     test_includeSimpleStartUpCost!()
+    test_includeBaseElasticDemand!()
+    test_includeVectorPrice!()
 end
 
 function _setup_common_variables() 
@@ -795,6 +797,40 @@ function test_includeSimpleStartUpCost!()
     register_tested_methods(includeSimpleStartUpCost!, 1)
 end
 
+function test_includeVectorPrice!()
+    # TODO: test for value::Dict
+    register_tested_methods(includeVectorPrice!, 1)
+end
+
+function test_includeBaseElasticDemand!()
+    # tests for value::Dict
+    (k, TL, LL) = _setup_common_variables()
+
+    d = Dict("Balance" => "PowerBalance_NO2", 
+             "Param" => "FirmDemand",
+             "price_elasticity" => 1.0, 
+             "normal_price" => 1.0 , 
+             "max_price" => 2.0, 
+             "min_price" => 1.0,
+             "threshold" => 0.15
+            )
+
+    bal = BaseBalance(
+            Id("a", "a"),
+            BaseCommodity(
+                Id("a", "a"), SequentialHorizon(1, Day(1))
+            )
+          )
+    par = MeanSeriesParam(ConstantTimeVector(1), ConstantTimeVector(1))
+    TL[Id("Balance", "PowerBalance_NO2")] = bal
+    LL[Id("Param", "FirmDemand")] = par
+
+    ret = includeBaseElasticDemand!(TL, LL, k, d)
+    _test_ret(ret, n=2, okvalue=true)
+    @test TL[Id(k.conceptname, k.instancename)] isa ElasticDemand
+    register_tested_methods(includeBaseElasticDemand!, 1)
+end
+
 function _test_ret(ret; n=0, okvalue=true, depstype=Vector{Id})
     @test ret isa Tuple{Bool, Any}
     (ok, deps) = ret