Allow expressions in Nodes state_initial and state_final parameters

assets/examples/56_final_state_decision_cyclic.iesopt.yaml

@@ -0,0 +1,73 @@
+parameters:
+    p_nom: 70
+    hours: 2.0
+    efficiency: 0.975
+    capex: 1000
+    opex_fom: 500
+    opex_vom: 1.8
+    T: 24
+
+config:
+  general:
+    version:
+      core: 2.6.2
+  optimization:
+    problem_type: LP
+    snapshots:
+      count: <T>
+    solver:
+      name: highs
+
+carriers:
+  electricity: {}
+
+components:
+  # --------------------------
+
+  storage:
+    type: Node
+    carrier: electricity
+    state_final: 0.5 * (<p_nom> + invest:value) * <hours>
+    state_cyclic: eq
+    has_state: true
+    state_lb: 0
+    state_ub: (<p_nom> + invest:value) * <hours>
+
+  node_electricity:
+    type: Node
+    carrier: electricity
+
+  # --------------------------
+
+  market_electricity:
+    type: Profile
+    carrier: electricity
+    mode: ranged
+    node_from: node_electricity
+    cost: [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]
+
+  # --------------------------
+
+  invest:
+    type: Decision
+    ub: 30
+    cost: ((<capex> + <opex_fom>) / 8760.0 * <T>)
+
+  # --------------------------
+
+  charge:
+    type: Unit
+    inputs: {electricity: node_electricity}
+    outputs: {electricity: storage}
+    conversion: 1 electricity -> <efficiency> electricity
+    capacity: (<p_nom> + invest:value) in:electricity
+    marginal_cost: <opex_vom> per out:electricity
+    objectives: {total_cost: 0.0 + (<p_nom> * <opex_fom>) / 8760.0 *<T>}
+
+  discharge:
+    type: Unit
+    inputs: {electricity: storage}
+    outputs: {electricity: node_electricity}
+    conversion: 1 electricity -> <efficiency> electricity
+    capacity: (<p_nom> + invest:value) out:electricity
+    marginal_cost: <opex_vom> per out:electricity

src/core/node.jl

@@ -65,15 +65,15 @@ balance equation. This allows using `Node`s for various storage tasks (like batt
     """
     state_cyclic::Symbol = :eq
 
-    raw"""```{"mandatory": "no", "values": "numeric", "unit": "energy", "default": "-"}```
-    Sets the initial state. Must be used in combination with `state_cyclic = disabled`.
+    raw"""```{"mandatory": "no", "values": "numeric, `decision:value`", "unit": "energy", "default": "-"}```
+    Sets the initial state. If both `state_initial` and `state_final` are set, `state_cyclic = disabled` is required.
     """
-    state_initial::_OptionalScalarInput = nothing
+    state_initial::Expression = @_default_expression(nothing)
 
-    raw"""```{"mandatory": "no", "values": "numeric", "unit": "energy", "default": "-"}```
-    Sets the final state. Must be used in combination with `state_cyclic = disabled`.
+    raw"""```{"mandatory": "no", "values": "numeric, `decision:value`", "unit": "energy", "default": "-"}```
+    Sets the final state. If both `state_initial` and `state_final` are set, `state_cyclic = disabled` is required.
     """
-    state_final::_OptionalScalarInput = nothing
+    state_final::Expression = @_default_expression(nothing)
 
     raw"""```{"mandatory": "no", "values": "``\\in [0, 1]``", "unit": "-", "default": "0"}```
     Per hour percentage loss of state (losing 1% should be set as `0.01`), will be scaled automatically for
@@ -128,10 +128,6 @@ end
 function _isvalid(node::Node)
     (node.state_cyclic in [:eq, :geq, :disabled]) || (@critical "<state_cyclic> invalid" node = node.name)
 
-    if !isnothing(node.state_final) && node.state_cyclic != :disabled
-        @critical "Nodes with a fixed final state need to set `state_cyclic` to `disabled`" node = node.name node.state_cyclic
-    end
-
     if !isnothing(node.etdf_group) && node.has_state
         @critical "Activating ETDF is not supported for stateful nodes" node = node.name
     end
@@ -210,6 +206,7 @@ include("node/var_state.jl")
 include("node/var_pf_theta.jl")
 include("node/con_state_bounds.jl")
 include("node/con_nodalbalance.jl")
+include("node/con_first_state.jl")
 include("node/con_last_state.jl")
 
 function _construct_expressions!(node::Node)
@@ -227,13 +224,29 @@ function _after_construct_variables!(node::Node)
     # We can now properly finalize the `state_lb`, and `state_ub`.
     _finalize(node.state_lb)
     _finalize(node.state_ub)
+    _finalize(node.state_initial)
+    _finalize(node.state_final)
+    # Check expressions
+    if node.state_initial.value isa Vector
+        @critical "[build] The initial value of a Node must be scalar." component = node.name
+    end
+    if node.state_final isa Vector
+        @critical "[build] The final value of a Node must be scalar." component = node.name
+    end
+    if node.state_cyclic !== :disabled && !isnothing(node.state_initial.value) && !isnothing(node.state_initial.value)
+        @critical(
+            "[build] `state_cyclic` has to be disabled if both `state_initial` and `state_final` are set.",
+            component = node.name
+        )
+    end
 
     return nothing
 end
 
 function _construct_constraints!(node::Node)
     _node_con_state_bounds!(node)
     _node_con_nodalbalance!(node) # 25% here
+    _node_con_first_state!(node)
     _node_con_last_state!(node)
     return nothing
 end

src/core/node/con_first_state.jl

@@ -0,0 +1,17 @@
+@doc raw"""
+    _node_con_last_state!(node::Node)
+
+Add the constraint the `node`'s state at the first Snapshot to the `model`, if
+`node.has_state == true` and `node.initial_state` is set.
+"""
+function _node_con_first_state!(node::Node)
+    node.has_state || return nothing
+    isnothing(node.state_initial.value) && return nothing
+    model = node.model
+    node.con.first_state = @constraint(
+        model,
+        node.var.state[1] == node.state_initial.value,
+        base_name = make_base_name(node, "first_state")
+    )
+    return nothing
+end

src/core/node/con_last_state.jl

@@ -45,16 +45,16 @@ function _node_con_last_state!(node::Node)
     lb = access(node.state_lb, t, OptionalScalarExpressionValue)
     ub = access(node.state_ub, t, OptionalScalarExpressionValue)
 
-    if !isnothing(node.state_final)
+    if !isnothing(node.state_final.value)
         # TODO: since lb/ub are potentially decision-containing expressions, we cannot check this here ...
         #       re-work this somehow (and consider that a user might want to force a final state that is out of "bounds", when using a decision)
         # if (!isnothing(lb) && (node.state_final < lb)) || (!isnothing(ub) && (node.state_final > ub))
         #     @warn "`state_final` is out of bounds and will be overwritten" node = node.name state_final =
         #         node.state_final lb ub
         # end
 
-        lb = node.state_final
-        ub = node.state_final
+        lb = node.state_final.value
+        ub = node.state_final.value
     end
 
     if !isnothing(lb) && node.nodal_balance != :create

src/core/node/var_state.jl

@@ -18,7 +18,5 @@ function _node_var_state!(node::Node)
 
     node.var.state = @variable(model, [t = get_T(model)], base_name = make_base_name(node, "state"), container = Array)
 
-    if !isnothing(node.state_initial)
-        JuMP.fix(node.var.state[1], node.state_initial; force=false)
-    end
+    return nothing
 end

src/parser.jl

@@ -292,6 +292,8 @@ function _parse_components!(model::JuMP.Model, @nospecialize(description::Dict{S
             # Convert to _Expression.
             state_lb = _convert_to_expression(model, pop!(prop, "state_lb", nothing))
             state_ub = _convert_to_expression(model, pop!(prop, "state_ub", nothing))
+            state_initial = _convert_to_expression(model, pop!(prop, "state_initial", nothing))
+            state_final = _convert_to_expression(model, pop!(prop, "state_final", nothing))
 
             # Convert to Symbol
             state_cyclic = Symbol(pop!(prop, "state_cyclic", :eq))
@@ -305,6 +307,8 @@ function _parse_components!(model::JuMP.Model, @nospecialize(description::Dict{S
                 soft_constraints_penalty,
                 state_lb,
                 state_ub,
+                state_initial,
+                state_final,
                 state_cyclic,
                 nodal_balance,
                 Dict(Symbol(k) => v for (k, v) in prop)...,

test/src/examples.jl

@@ -276,3 +276,10 @@ end
     @test JuMP.value.(get_component(model, "pipeline.invest").var.value) ≈ 42.0 atol = 1e-2
     @test get_component(model, "pipeline.invest").cost ≈ 60067.3621 atol = 1e-2
 end
+
+@testitem "56_final_state_decision_cyclic" tags = [:examples] setup = [Dependencies, TestExampleModule] begin
+    model = TestExampleModule.check(; obj=-1.7867962768e+03)
+
+    @test JuMP.value(get_component(model, "storage").var.state[1]) == 100
+    @test sum(JuMP.value, get_component(model, "storage").exp.injection) == 0
+end