Updates for IRTOC runs on EI and update to PSI 0.34

.gitignore

@@ -121,5 +121,8 @@ $RECYCLE.BIN/
 # Windows shortcuts
 *.lnk
 
+# Claude - modify to include claude files 
+.claude/
+
 ## Acknowledgements
 # Many thanks to `https://gitignore.io/`, written and maintained by Joe Blau, which contributed much material to this gitignore file.

Project.toml

@@ -16,14 +16,14 @@ PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 
 [compat]
-DataStructures = "^0.18"
+DataStructures = "~0.18, ~0.19"
 Dates = "1"
 DocStringExtensions = "~0.8, ~0.9"
-InfrastructureSystems = "2"
+InfrastructureSystems = "3"
 JuMP = "1"
 MPI = "^0.20"
 MathOptInterface = "1"
-PowerNetworkMatrices = "0.12"
-PowerSystems = "4"
-PowerSimulations = "^0.30"
+PowerNetworkMatrices = "^0.20"
+PowerSystems = "5.3"
+PowerSimulations = "0.34"
 julia = "^1.6"

src/algorithms/sequential_algorithm.jl

@@ -20,21 +20,14 @@ function build_main_problem!(
     template::MultiProblemTemplate,
     sys::PSY.System,
 )
-    branch_models_dict = keys(PSI.get_branch_models(template))
     device_models_dict = keys(PSI.get_device_models(template))
-    if :TwoTerminalHVDCLine ∈ branch_models_dict ||
-       :TwoTerminalVSCDCLine ∈ branch_models_dict
-        has_hvdc = true
-    else
-        has_hvdc = false
-    end
     for k in keys(container.subproblems)
-        subsystem_buses = PSY.get_components(PSY.ACBus, sys; subsystem_name=k)
+        subsystem_buses = PSY.get_available_components(PSY.ACBus, sys; subsystem_name=k)
         subsystem_bus_nos = [PSY.get_number(b) for b in subsystem_buses]
         container.subproblem_bus_map[k] = subsystem_bus_nos
         subsystem_static_injectors =
-            PSY.get_components(PSY.StaticInjection, sys; subsystem_name=k)
-        subsystem_buses = PSY.get_components(PSY.ACBus, sys; subsystem_name=k)
+            PSY.get_available_components(PSY.StaticInjection, sys; subsystem_name=k)
+        subsystem_buses = PSY.get_available_components(PSY.ACBus, sys; subsystem_name=k)
         for si in subsystem_static_injectors
             if Symbol(typeof(si)) ∈ device_models_dict
                 if !(PSY.get_bus(si) ∈ subsystem_buses)
@@ -46,12 +39,9 @@ function build_main_problem!(
                 end
             end
         end
-        subsystem_hvdcs = PSY.get_components(
-            Union{PSY.TwoTerminalHVDCLine, PSY.TwoTerminalVSCDCLine},
-            sys;
-            subsystem_name=k,
-        )
-        if has_hvdc
+        subsystem_hvdcs =
+            PSY.get_available_components(PSY.TwoTerminalHVDC, sys; subsystem_name=k)
+        if _has_hvdc_model(template)
             for hvdc in subsystem_hvdcs
                 from_bus_no = PSY.get_number(PSY.get_from(PSY.get_arc(hvdc)))
                 to_bus_no = PSY.get_number(PSY.get_to(PSY.get_arc(hvdc)))
@@ -67,6 +57,20 @@ function build_main_problem!(
     end
 end
 
+function _has_hvdc_model(template::MultiProblemTemplate)
+    branch_models_dict = keys(PSI.get_branch_models(template))
+    for hvdc_type in CONCRETE_HVDC_TYPES
+        if hvdc_type in branch_models_dict
+            return true
+        end
+    end
+    return false
+end
+
+# Note: With the addition of 3D results processing, we can eliminate this design of writing
+# all subsystem results to a signle container. This way we avoid overwriting the results from 
+# different subsystem problems and can still access all relevant results using the API. 
+
 # The drawback of this approach is that it will loop over the results twice
 # once to write into the main container and a second time when writing into the
 # store. The upside of this approach is that doesn't require overloading write_model_XXX_results!
@@ -79,23 +83,30 @@ function write_results_to_main_container(container::MultiOptimizationContainer)
             subproblem_data_field = getproperty(subproblem, field)
             main_container_data_field = getproperty(container, field)
             for (key, src) in subproblem_data_field
+                @debug "writing $key to main container"
                 if src isa JuMP.Containers.SparseAxisArray
                     @debug "Skip SparseAxisArray" field key
                     continue
                 end
                 num_dims = ndims(src)
-                num_dims > 2 && error("ndims = $(num_dims) is not supported yet")
-                data = nothing
-                data = PSI.jump_value.(src)
                 dst = main_container_data_field[key]
                 if num_dims == 1
-                    dst[1:length(axes(src)[1])] = data
+                    columns = _get_main_container_columns(container, k, key, src)
+                    for col in columns
+                        if isassigned(src, col)
+                            dst[col] = PSI.jump_value.(src[col])
+                        end
+                    end
                 elseif num_dims == 2
                     columns = _get_main_container_columns(container, k, key, src)
-                    len = length(axes(src)[2])
-                    dst[columns, 1:len] = PSI.jump_value.(src[columns, :])
+                    for col in columns
+                        for t in axes(src)[2]
+                            if isassigned(src, col, t)
+                                dst[col, t] = PSI.jump_value.(src[col, t])
+                            end
+                        end
+                    end
                 elseif num_dims == 3
-                    # TODO: untested
                     axis1 = axes(src)[1]
                     axis2 = axes(src)[2]
                     len = length(axes(src)[3])

src/core/definitions.jl

@@ -1,3 +1,6 @@
 const CONTAINER_FIELDS = [:variables, :aux_variables, :constraints, :expressions, :duals]
 const ALL_CONTAINER_FIELDS =
     [:variables, :aux_variables, :constraints, :expressions, :duals, :parameters]
+
+const CONCRETE_HVDC_TYPES =
+    [:TwoTerminalGenericHVDCLine, :TwoTerminalLCCLine, :TwoTerminalGenericHVDCLine]

src/models/branch_models.jl

@@ -6,13 +6,6 @@ function PSI.construct_device!(
     network_model::PSI.NetworkModel{<:PSI.AbstractPTDFModel},
 )
     devices = PSI.get_available_components(model, sys)
-    PSI.add_variables!(
-        container,
-        PSI.FlowActivePowerVariable,
-        network_model,
-        devices,
-        PSI.StaticBranchUnbounded(),
-    )
     PSI.add_parameters!(container, StateEstimationFlows, devices, model)
     PSI.add_feedforward_arguments!(container, model, devices)
     return
@@ -30,7 +23,7 @@ function PSI.add_parameters!(
         container,
         StateEstimationFlows(),
         D,
-        ISOPT.VariableKey{PSI.FlowActivePowerVariable, D}(""),
+        ISOPT.ExpressionKey{PSI.PTDFBranchFlow, D}(""),
         branch_names,
         time_steps,
     )
@@ -54,7 +47,7 @@ function PSI._make_flow_expressions!(
 )
     # NOTE - can alternatively avoid passing the system by ensuring the parameter container for state estimation injections
     # is ordered by bus number as is the case for the Active Power Balance expressions
-    all_buses = PSY.get_components(
+    all_buses = PSY.get_available_components(
         x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
         PSY.ACBus,
         sys;
@@ -118,20 +111,14 @@ function PSI._make_flow_expressions!(
     container::PSI.OptimizationContainer,
     branches::Vector{String},
     time_steps::UnitRange{Int},
-    ptdf::PSI.ValidPTDFS,
-    nodal_balance_expressions::PSI.JuMPAffineExpressionDArray,
+    ptdf::Union{PNM.PTDF, PNM.VirtualPTDF},
+    nodal_balance_expressions::PSI.JuMPAffineExpressionDArrayIntInt,
     state_estimation_injections,
     state_estimation_flows,
     branch_Type::DataType,
     sys::PSY.System,
 )
-    branch_flow_expr = PSI.add_expression_container!(
-        container,
-        PSI.PTDFBranchFlow(),
-        branch_Type,
-        branches,
-        time_steps,
-    )
+    branch_flow_expr = PSI.get_expression(container, PSI.PTDFBranchFlow(), branch_Type)
 
     jump_model = PSI.get_jump_model(container)
 
@@ -183,8 +170,8 @@ function PSI.add_constraints!(
     ptdf = PSI.get_PTDF_matrix(network_model)
     # This is a workaround to not call the same list comprehension to find
     # The subset of branches of type B in the PTDF
-    flow_variables = PSI.get_variable(container, PSI.FlowActivePowerVariable(), B)
-    branches = flow_variables.axes[1]
+    flows = PSI.get_expression(container, PSI.PTDFBranchFlow(), B)
+    branches = flows.axes[1]
     time_steps = PSI.get_time_steps(container)
     branch_flow = PSI.add_constraints_container!(
         container,
@@ -199,7 +186,6 @@ function PSI.add_constraints!(
         PSI.get_parameter_array(container, StateEstimationInjections(), PSY.ACBus)
 
     state_estimation_flows = PSI.get_parameter(container, StateEstimationFlows(), B)
-    flow_variables = PSI.get_variable(container, PSI.FlowActivePowerVariable(), B)
     branch_flow_expr = PSI._make_flow_expressions!(
         container,
         branches,
@@ -216,7 +202,7 @@ function PSI.add_constraints!(
         for t in time_steps
             branch_flow[name, t] = JuMP.@constraint(
                 jump_model,
-                branch_flow_expr[name, t] - flow_variables[name, t] == 0.0
+                branch_flow_expr[name, t] - flows[name, t] == 0.0
             )
         end
     end
@@ -231,6 +217,7 @@ function PSI.construct_device!(
     network_model::PSI.NetworkModel{<:PSI.AbstractPTDFModel},
 )
     devices = PSI.get_available_components(model, sys)
+    PSI.add_expressions!(container, PSI.PTDFBranchFlow, devices, model, network_model)
     # NOTE - changes required in handling of feedforwards? 
     PSI.add_feedforward_constraints!(container, model, devices)
     PSI.add_constraints!(

src/models/network_models.jl

@@ -42,7 +42,6 @@ function PSI.construct_network!(
         model,
     )
     PSI.add_constraints!(container, PSI.CopperPlateBalanceConstraint, sys, model)
-    PSI.add_constraints!(container, PSI.NodalBalanceActiveConstraint, sys, model)
     PSI.add_constraint_dual!(container, sys, model)
     return
 end
@@ -79,7 +78,7 @@ function PSI.add_to_expression!(
     parameter_array =
         PSI.get_parameter_array(container, StateEstimationInjections(), PSY.ACBus)
     subsys = PSI.get_subsystem(network_model)
-    all_buses = PSY.get_components(
+    all_buses = PSY.get_available_components(
         x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
         PSY.ACBus,
         sys;
@@ -88,12 +87,12 @@ function PSI.add_to_expression!(
     # These are the buses not in the same subsystem as the one being built
 
     expression = PSI.get_expression(container, PSI.ActivePowerBalance(), PSY.ACBus)
-    radial_network_reduction = PSI.get_radial_network_reduction(network_model)
+    nrd = PSI.get_network_reduction(network_model)
     for b in all_buses, t in PSI.get_time_steps(container)
         if PSY.has_component(sys, subsys, b)
             continue
         end
-        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, b)
+        bus_no = PNM.get_mapped_bus_number(nrd, b)
         PSI._add_to_jump_expression!(
             expression[bus_no, t],
             parameter_array[string(bus_no), t],
@@ -130,13 +129,7 @@ function PSI.add_parameters!(
     time_steps = PSI.get_time_steps(container)
     subsys = PSI.get_subsystem(network_model)
 
-    all_buses = PSY.get_components(
-        x -> PSY.get_bustype(x) != PSY.ACBusTypes.ISOLATED,
-        PSY.ACBus,
-        sys;
-    )
-
-    bus_numbers = [string(PSY.get_number(b)) for b in all_buses]
+    bus_numbers = string.(PNM.get_bus_axis(PSI.get_PTDF_matrix(network_model)))
     @assert !isempty(bus_numbers)
 
     parameter_container = PSI.add_param_container!(
@@ -162,6 +155,7 @@ function PSI.initialize_system_expressions!(
     container::PSI.OptimizationContainer,
     network_model::PSI.NetworkModel{SplitAreaPTDFPowerModel},
     subnetworks::Dict{Int, Set{Int}},
+    ::PSI.BranchModelContainer,
     system::PSY.System,
     bus_reduction_map::Dict{Int64, Set{Int64}},
 )
@@ -198,12 +192,12 @@ function PSI.add_to_expression!(
     variable = PSI.get_variable(container, U(), V)
     area_expr = PSI.get_expression(container, T(), PSY.Area)
     nodal_expr = PSI.get_expression(container, T(), PSY.ACBus)
-    radial_network_reduction = PSI.get_radial_network_reduction(network_model)
+    nrd = PSI.get_network_reduction(network_model)
     for d in devices
         name = PSY.get_name(d)
         device_bus = PSY.get_bus(d)
         area_name = PSY.get_name(PSY.get_area(device_bus))
-        bus_no = PNM.get_mapped_bus_number(radial_network_reduction, device_bus)
+        bus_no = PNM.get_mapped_bus_number(nrd, device_bus)
         for t in PSI.get_time_steps(container)
             PSI._add_to_jump_expression!(
                 area_expr[area_name, t],
@@ -281,19 +275,17 @@ function _update_parameter_values!(
 ) where {T <: Union{JuMP.VariableRef, Float64}}
     state = PSI.get_system_states(simulation_state)
     state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
-
     if !isfinite(first(state_values))
         @warn "first value not present, updating state estimation injections from decision state"
         state = PSI.get_decision_states(simulation_state)
         state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
+    elseif size(parameter_array)[2] > size(state_values)[2]
+        @warn "Cannot update; state estimation injection parameter has more timesteps than the system state, updating state estimation injections from decision state"
+        state = PSI.get_decision_states(simulation_state)
+        state_values = PSI.get_dataset_values(state, PSI.get_attribute_key(attributes))
+    else
+        @info "Updating state estimation injection parameters from system state"
     end
-
-    if size(parameter_array)[2] > size(state_values)[2]
-        error(
-            "Cannot update: state estimation injection parameter has more timesteps than the state used for updating.",
-        )
-    end
-
     component_names, time = axes(parameter_array)
     for t in time
         for name in component_names

src/multi_optimization_container.jl

@@ -6,7 +6,6 @@ Base.@kwdef mutable struct MultiOptimizationContainer{T <: DecompositionAlgorith
     time_steps::UnitRange{Int}
     resolution::Dates.TimePeriod
     settings::PSI.Settings
-    settings_copy::PSI.Settings
     variables::Dict{ISOPT.VariableKey, AbstractArray}
     aux_variables::Dict{ISOPT.AuxVarKey, AbstractArray}
     duals::Dict{ISOPT.ConstraintKey, AbstractArray}
@@ -59,7 +58,6 @@ function MultiOptimizationContainer(
         time_steps=1:1,
         resolution=IS.time_period_conversion(resolution),
         settings=settings,
-        settings_copy=PSI.copy_for_serialization(settings),
         variables=Dict{ISOPT.VariableKey, AbstractArray}(),
         aux_variables=Dict{ISOPT.AuxVarKey, AbstractArray}(),
         duals=Dict{ISOPT.ConstraintKey, AbstractArray}(),
@@ -187,3 +185,18 @@ function PSI.serialize_optimization_model(
     container::MultiOptimizationContainer,
     save_path::String,
 ) end
+
+function PSI.get_column_names(
+    ::MultiOptimizationContainer,
+    field::Symbol,
+    subcontainer,
+    key::PSI.OptimizationContainerKey,
+)
+    return if field == :parameters
+        # Parameters are stored in ParameterContainer.
+        PSI.get_column_names(key, subcontainer)
+    else
+        # The others are in DenseAxisArrays.
+        PSI.get_column_names_from_axis_array(key, subcontainer)
+    end
+end