From dfaefddc7b38e9158411823dece0430cb42337d7 Mon Sep 17 00:00:00 2001
From: Greg Schivley <greg.schivley@gmail.com>
Date: Tue, 12 May 2026 10:01:17 -0400
Subject: [PATCH] Enhance output writing functions with caching and
 transposition improvements

- Added optional `cache` argument to multiple output writing functions (`write_curtailment`, `write_emissions`, `write_energy_revenue`, `write_power`, `write_storage`, `write_nse`, `write_price`, `write_storagedual`, etc.) to allow reuse of extracted model outputs, reducing memory allocations.
- Implemented `resource_time_scratch!` and `scaled_resource_time_matrix!` functions to handle resource time matrices efficiently with caching.
- Updated CSV writing functions to use `write_transposed_csv` for improved performance and consistency.
- Introduced a new test suite for output caching and transposition, ensuring correctness and performance of the new features.
- Refactored existing code to utilize the new caching mechanism, enhancing overall efficiency.

Co-authored-by: Copilot <copilot@github.com>
---
 CHANGELOG.md                                  |   1 +
 .../full_time_series_reconstruction.jl        |  35 ++-
 .../write_reserve_margin_slack.jl             |   4 +-
 src/write_outputs/dftranspose.jl              |  51 ++--
 .../write_hourly_matching_prices.jl           |   5 +-
 .../write_hourly_matching_slack.jl            |   4 +-
 .../write_opwrap_lds_dstor.jl                 |   4 +-
 .../write_opwrap_lds_stor_init.jl             |   6 +-
 src/write_outputs/output_cache.jl             | 181 ++++++++++++
 src/write_outputs/reserves/write_rsv.jl       |   4 +-
 .../transmission/write_transmission_flows.jl  |   2 +-
 .../transmission/write_transmission_losses.jl |   4 +-
 src/write_outputs/write_angles.jl             |   4 +-
 src/write_outputs/write_capacityfactor.jl     |  17 +-
 src/write_outputs/write_charge.jl             |  34 ++-
 src/write_outputs/write_charging_cost.jl      |  37 ++-
 src/write_outputs/write_curtailment.jl        |  23 +-
 src/write_outputs/write_emissions.jl          |  32 ++-
 src/write_outputs/write_energy_revenue.jl     |  29 +-
 src/write_outputs/write_fuel_consumption.jl   |   5 +-
 src/write_outputs/write_nse.jl                |  19 +-
 src/write_outputs/write_outputs.jl            |  64 +++--
 src/write_outputs/write_power.jl              |  17 +-
 src/write_outputs/write_power_balance.jl      |   4 +-
 src/write_outputs/write_price.jl              |  18 +-
 src/write_outputs/write_reliability.jl        |   4 +-
 src/write_outputs/write_storage.jl            |  36 ++-
 src/write_outputs/write_storagedual.jl        |  33 ++-
 test/runtests.jl                              |   4 +
 test/test_output_cache.jl                     | 260 ++++++++++++++++++
 30 files changed, 750 insertions(+), 191 deletions(-)
 create mode 100644 src/write_outputs/output_cache.jl
 create mode 100644 test/test_output_cache.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e32e44fe5d..7f149a4f1a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 - Fix writing of net revenue to include all sources of revenue, not just energy revenue (#855).
+- Reduce output-writing memory spikes by caching extracted model values, reusing large scratch arrays, and avoiding copy-heavy transposed DataFrame writes (#908).
 
 ## [0.4.6] - 2026-01-06
 
diff --git a/src/time_domain_reduction/full_time_series_reconstruction.jl b/src/time_domain_reduction/full_time_series_reconstruction.jl
index 3b7900bed5..423bad569e 100644
--- a/src/time_domain_reduction/full_time_series_reconstruction.jl
+++ b/src/time_domain_reduction/full_time_series_reconstruction.jl
@@ -38,34 +38,39 @@ function full_time_series_reconstruction(
 
     # Get a matrix of the input DataFrame
     DFMatrix = Matrix(DF)
-    # Initialize an array to add the reconstructed data to
-    recon = ["t$t" for t in 1:(TimestepsPerRepPeriod * numPeriods)]
+    total_timesteps = TimestepsPerRepPeriod * numPeriods
 
     # Find the index of the row with the first time step
     t1 = findfirst(x -> x == "t1", DF[!, 1])
 
     # Reconstruction of all hours of the year from TDR
-    for j in range(2, ncol(DF))
-        col = DF[t1:end, j]
-        recon_col = []
-        for i in range(1, numPeriods)
+    recon = Matrix{Any}(undef, total_timesteps, ncol(DF))
+    recon[:, 1] = ["t$t" for t in 1:total_timesteps]
+    for j in 2:ncol(DF)
+        col = DFMatrix[t1:end, j]
+        next_row = 1
+        for i in 1:numPeriods
             index = Period_map[i, "Rep_Period_Index"]
-            recon_temp = col[(TimestepsPerRepPeriod * index - (TimestepsPerRepPeriod - 1)):(TimestepsPerRepPeriod * index)]
-            recon_col = [recon_col; recon_temp]
+            source_start = TimestepsPerRepPeriod * (index - 1) + 1
+            source_end = TimestepsPerRepPeriod * index
+            target_end = next_row + TimestepsPerRepPeriod - 1
+            recon[next_row:target_end, j] = col[source_start:source_end]
+            next_row = target_end + 1
         end
-        recon = [recon recon_col]
     end
-    reconDF = DataFrame(recon, :auto)
+    reconDF = DataFrame(recon, names(DF))
 
     # Insert rows that were above "t1" in the original DataFrame (e.g. "Zone" and "AnnualSum") if present
-    for i in range(1, t1 - 1)
-        insert!(reconDF, i, DFMatrix[i, 1:end], promote = true)
+    if t1 > 1
+        reconDF = vcat(DataFrame(DFMatrix[1:(t1 - 1), :], names(DF)), reconDF)
     end
 
     # Repeat the last rows of the year to fill in the gap (should be 24 hours for non-leap year)
     end_diff = WeightTotal - nrow(reconDF) + 1
-    new_rows = reconDF[(nrow(reconDF) - end_diff):nrow(reconDF), 1:end]
-    new_rows[!, 1] = ["t$t" for t in (WeightTotal - end_diff):WeightTotal]
-    reconDF = [reconDF; new_rows]
+    if end_diff > 0
+        new_rows = copy(reconDF[(nrow(reconDF) - end_diff):nrow(reconDF), 1:end])
+        new_rows[!, 1] = ["t$t" for t in (WeightTotal - end_diff):WeightTotal]
+        append!(reconDF, new_rows)
+    end
     return reconDF
 end
diff --git a/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl b/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
index 3248353b65..d29ac6d84c 100644
--- a/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
+++ b/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
@@ -22,8 +22,8 @@ function write_reserve_margin_slack(path::AbstractString,
         end
         dfResMar_slack = hcat(dfResMar_slack,
             DataFrame(temp_ResMar_slack, [Symbol("t$t") for t in 1:T]))
-        CSV.write(joinpath(path, "ReserveMargin_slack_and_penalties.csv"),
-            dftranspose(dfResMar_slack, false),
+        write_transposed_csv(joinpath(path, "ReserveMargin_slack_and_penalties.csv"),
+            dfResMar_slack,
             writeheader = false)
     end
     return nothing
diff --git a/src/write_outputs/dftranspose.jl b/src/write_outputs/dftranspose.jl
index ec9d1a5f39..a923a2d79d 100644
--- a/src/write_outputs/dftranspose.jl
+++ b/src/write_outputs/dftranspose.jl
@@ -1,26 +1,37 @@
 ################################################################################
-## function dftranspose(df)
-##
-## inputs: df - [DataFrame] a DataFrame object to be transposed
-## results t - [DataFrame] a transposed version of the df DataFrame.
-##		   withhead - [Boolean] if True, first column of df will become column
-##		   names for t. Otherwise, first column will first row and column names
-##		   will be generic (e.g. x1:xN)
-##
-## Note this function is necessary because no stock function to transpose
-## DataFrames appears to exist.
+## helpers for writing transposed output tables
 ################################################################################
 @doc raw"""
-   df = dftranspose(df::DataFrame, withhead::Bool)
+    transpose_output_dataframe(df::DataFrame; withhead::Bool=false)
 
-Returns a transpose of a Dataframe.
+Return a transposed copy of an output DataFrame without collecting each row into
+an intermediate vector.
 """
-function dftranspose(df::DataFrame, withhead::Bool)
-    if withhead
-        colnames = cat(:Row, Symbol.(df[!, 1]), dims = 1)
-        return DataFrame([[names(df)]; collect.(eachrow(df))], colnames)
-    else
-        return DataFrame([[names(df)]; collect.(eachrow(df))],
-            [:Row; Symbol.("x", axes(df, 1))])
+function transpose_output_dataframe(df::DataFrame; withhead::Bool = false)
+    row_count, col_count = size(df)
+    colnames = withhead ? Symbol[:Row; Symbol.(df[!, 1])] :
+               Symbol[:Row; Symbol.("x", 1:row_count)]
+    transposed = Matrix{Any}(undef, col_count, row_count + 1)
+    headers = names(df)
+
+    for col in 1:col_count
+        transposed[col, 1] = headers[col]
+    end
+
+    for row in 1:row_count
+        for col in 1:col_count
+            transposed[col, row + 1] = df[row, col]
+        end
+    end
+
+    return DataFrame(transposed, colnames)
+end
+
+function write_transposed_csv(filename::AbstractString, df::DataFrame; kwargs...)
+    options = Dict{Symbol, Any}(pairs(kwargs))
+    if haskey(options, :writeheader) && !haskey(options, :header)
+        options[:header] = options[:writeheader]
+        delete!(options, :writeheader)
     end
-end # End dftranpose()
+    return CSV.write(filename, transpose_output_dataframe(df); options...)
+end
diff --git a/src/write_outputs/hourly_matching_requirement/write_hourly_matching_prices.jl b/src/write_outputs/hourly_matching_requirement/write_hourly_matching_prices.jl
index 7a6d209564..30330d2d11 100644
--- a/src/write_outputs/hourly_matching_requirement/write_hourly_matching_prices.jl
+++ b/src/write_outputs/hourly_matching_requirement/write_hourly_matching_prices.jl
@@ -23,7 +23,8 @@ function write_hourly_matching_prices(path::AbstractString,
     auxNew_Names = [Symbol("Zone"); [Symbol("t$t") for t in 1:T]]
     rename!(dfHourlyMatchPrices, auxNew_Names)
 
-    CSV.write(joinpath(path, "hourly_matching_prices.csv"),
-        dftranspose(dfHourlyMatchPrices, false), writeheader = false)
+    write_transposed_csv(joinpath(path, "hourly_matching_prices.csv"),
+        dfHourlyMatchPrices,
+        writeheader = false)
     return nothing
 end
diff --git a/src/write_outputs/hourly_matching_requirement/write_hourly_matching_slack.jl b/src/write_outputs/hourly_matching_requirement/write_hourly_matching_slack.jl
index 7cc14ab54c..968981b636 100644
--- a/src/write_outputs/hourly_matching_requirement/write_hourly_matching_slack.jl
+++ b/src/write_outputs/hourly_matching_requirement/write_hourly_matching_slack.jl
@@ -28,8 +28,8 @@ function write_hourly_matching_slack(path::AbstractString,
         end
         dfHM_slack = hcat(dfHM_slack,
             DataFrame(temp_HM_slack, [Symbol("t$t") for t in 1:T]))
-        CSV.write(joinpath(path, "HourlyMatching_slack_and_penalties.csv"),
-            dftranspose(dfHM_slack, false),
+        write_transposed_csv(joinpath(path, "HourlyMatching_slack_and_penalties.csv"),
+            dfHM_slack,
             writeheader = false)
     end
     return nothing
diff --git a/src/write_outputs/long_duration_storage/write_opwrap_lds_dstor.jl b/src/write_outputs/long_duration_storage/write_opwrap_lds_dstor.jl
index 875d8e6f86..0b943724c6 100644
--- a/src/write_outputs/long_duration_storage/write_opwrap_lds_dstor.jl
+++ b/src/write_outputs/long_duration_storage/write_opwrap_lds_dstor.jl
@@ -26,7 +26,5 @@ function write_opwrap_lds_dstor(path::AbstractString, inputs::Dict, setup::Dict,
     dfdStorage = hcat(dfdStorage, DataFrame(dsoc, :auto))
     auxNew_Names = [Symbol("Resource"); Symbol("Zone"); [Symbol("w$t") for t in 1:W]]
     rename!(dfdStorage, auxNew_Names)
-    CSV.write(joinpath(path, "dStorage.csv"),
-        dftranspose(dfdStorage, false),
-        header = false)
+    write_transposed_csv(joinpath(path, "dStorage.csv"), dfdStorage, header = false)
 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
index e0a731d177..3455bfd4aa 100644
--- 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
@@ -29,8 +29,8 @@ function write_opwrap_lds_stor_init(path::AbstractString,
     dfStorageInit = hcat(dfStorageInit, DataFrame(socw, :auto))
     auxNew_Names = [Symbol("Resource"); Symbol("Zone"); [Symbol("n$t") for t in 1:NPeriods]]
     rename!(dfStorageInit, auxNew_Names)
-    CSV.write(joinpath(path, "StorageInit.csv"),
-        dftranspose(dfStorageInit, false),
+    write_transposed_csv(joinpath(path, "StorageInit.csv"),
+        dfStorageInit,
         header = false)
 
     # Write storage evolution over full time horizon
@@ -80,6 +80,6 @@ function write_opwrap_lds_stor_init(path::AbstractString,
     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)
+    write_transposed_csv(joinpath(path, "StorageEvol.csv"), df_SOC_t, writeheader = false)
 
 end
diff --git a/src/write_outputs/output_cache.jl b/src/write_outputs/output_cache.jl
new file mode 100644
index 0000000000..13f7e829d0
--- /dev/null
+++ b/src/write_outputs/output_cache.jl
@@ -0,0 +1,181 @@
+@doc raw"""
+        OutputCache
+
+Container for extracted model outputs reused across write-output functions to
+reduce repeated JuMP value materialization and large temporary allocations.
+
+# Fields
+- `scale_factor::Float64`: Output scaling multiplier (1 or `ModelScalingFactor`).
+- `resource_time_scratch::Matrix{Float64}`: Reusable `G x T` workspace.
+- `price::Union{Nothing, Matrix{Float64}}`: Cached locational marginal prices (`T x Z`).
+- `vP::Matrix{Float64}`: Cached dispatch (`G x T`).
+- `eTotalCap::Vector{Float64}`: Cached endogenously available capacity (`G`).
+- `vCHARGE`, `vCHARGE_FLEX`, `vUSE`, `vCHARGE_VRE_STOR`, `vCHARGE_ALLAM`:
+    optional cached charging/consumption matrices.
+- `vS`, `vS_HYDRO`, `vS_FLEX`, `vS_VRE_STOR`: optional cached storage-state matrices.
+- `vNSE::Union{Nothing, Array{Float64,3}}`: optional cached non-served energy tensor.
+- `eEmissionsByZone::Union{Nothing, Matrix{Float64}}`: optional cached zonal emissions.
+"""
+struct OutputCache
+    scale_factor::Float64
+    resource_time_scratch::Matrix{Float64}
+    price::Union{Nothing, Matrix{Float64}}
+    vP::Matrix{Float64}
+    eTotalCap::Vector{Float64}
+    vCHARGE::Union{Nothing, Matrix{Float64}}
+    vCHARGE_FLEX::Union{Nothing, Matrix{Float64}}
+    vUSE::Union{Nothing, Matrix{Float64}}
+    vCHARGE_VRE_STOR::Union{Nothing, Matrix{Float64}}
+    vCHARGE_ALLAM::Union{Nothing, Matrix{Float64}}
+    vS::Union{Nothing, Matrix{Float64}}
+    vS_HYDRO::Union{Nothing, Matrix{Float64}}
+    vS_FLEX::Union{Nothing, Matrix{Float64}}
+    vS_VRE_STOR::Union{Nothing, Matrix{Float64}}
+    vNSE::Union{Nothing, Array{Float64, 3}}
+    eEmissionsByZone::Union{Nothing, Matrix{Float64}}
+end
+
+_extract_output_data(var::JuMP.Containers.DenseAxisArray) = var.data
+_extract_output_data(var::AbstractArray) = var
+
+_extract_output_matrix(var)::Matrix{Float64} = Matrix{Float64}(Array(value.(
+    _extract_output_data(var))))
+
+_extract_output_vector(var)::Vector{Float64} = vec(Array(value.(
+    _extract_output_data(var))))
+
+@doc raw"""
+    build_output_cache(
+        EP::Model,
+        inputs::Dict,
+        setup::Dict,
+        output_settings_d::Dict = setup["WriteOutputsSettingsDict"];
+        selective::Bool = false)
+
+Build an `OutputCache` for write-output routines.
+
+When `selective=true`, only values needed by enabled outputs in
+`output_settings_d` are materialized; otherwise, common output matrices are
+eagerly extracted once.
+
+# Arguments
+- `EP::Model`: Solved JuMP model.
+- `inputs::Dict`: Parsed GenX input data.
+- `setup::Dict`: Run settings.
+- `output_settings_d::Dict`: Output toggles.
+- `selective::Bool=false`: Enable selective extraction.
+
+# Returns
+- `OutputCache`: Cache object shared by write-output functions.
+"""
+function build_output_cache(
+    EP::Model,
+    inputs::Dict,
+    setup::Dict,
+    output_settings_d::Dict = setup["WriteOutputsSettingsDict"];
+    selective::Bool = false)
+    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1.0
+    G = inputs["G"]
+    T = inputs["T"]
+
+    needs_power = output_settings_d["WritePower"] ||
+                  output_settings_d["WriteCapacityFactor"] ||
+                  output_settings_d["WriteCurtailment"] ||
+                  output_settings_d["WriteEnergyRevenue"] ||
+                  output_settings_d["WriteChargingCost"] ||
+                  output_settings_d["WriteNetRevenue"]
+    needs_capacity = output_settings_d["WriteCapacityFactor"] ||
+                     output_settings_d["WriteCurtailment"]
+    needs_charge = output_settings_d["WriteCharge"] ||
+                   output_settings_d["WriteChargingCost"] ||
+                   output_settings_d["WriteNetRevenue"]
+    needs_storage = output_settings_d["WriteStorage"] ||
+                    output_settings_d["WriteStorageDual"] ||
+                    output_settings_d["WriteNetRevenue"]
+    needs_price = has_duals(EP) == 1 && (
+        output_settings_d["WritePrice"] ||
+        output_settings_d["WriteEnergyRevenue"] ||
+        output_settings_d["WriteChargingCost"] ||
+        output_settings_d["WriteNetRevenue"])
+    needs_scratch = needs_power || needs_capacity || needs_storage ||
+                    output_settings_d["WriteNSE"] ||
+                    output_settings_d["WriteEmissions"] ||
+                    output_settings_d["WriteEnergyRevenue"] ||
+                    output_settings_d["WriteChargingCost"]
+
+    return OutputCache(
+        scale_factor,
+        selective && !needs_scratch ? zeros(Float64, 0, 0) : zeros(Float64, G, T),
+        needs_price ? locational_marginal_price(EP, inputs, setup) : nothing,
+        selective && !needs_power ? zeros(Float64, 0, 0) : _extract_output_matrix(EP[:vP]),
+        selective && !needs_capacity ? zeros(Float64, 0) : _extract_output_vector(EP[:eTotalCap]),
+        selective && !needs_charge || isempty(inputs["STOR_ALL"]) ? nothing : _extract_output_matrix(EP[:vCHARGE]),
+        selective && !needs_charge || isempty(inputs["FLEX"]) ? nothing : _extract_output_matrix(EP[:vCHARGE_FLEX]),
+        selective && !needs_charge || isempty(inputs["ELECTROLYZER"]) ? nothing : _extract_output_matrix(EP[:vUSE]),
+        selective && !needs_charge || isempty(inputs["VRE_STOR"]) ? nothing : _extract_output_matrix(EP[:vCHARGE_VRE_STOR]),
+        selective && !needs_charge || isempty(inputs["ALLAM_CYCLE_LOX"]) ? nothing : _extract_output_matrix(EP[:vCHARGE_ALLAM]),
+        selective && !needs_storage || isempty(inputs["STOR_ALL"]) ? nothing : _extract_output_matrix(EP[:vS]),
+        selective && !needs_storage || isempty(inputs["HYDRO_RES"]) ? nothing : _extract_output_matrix(EP[:vS_HYDRO]),
+        selective && !needs_storage || isempty(inputs["FLEX"]) ? nothing : _extract_output_matrix(EP[:vS_FLEX]),
+        selective && !needs_storage || isempty(inputs["VRE_STOR"]) ? nothing : _extract_output_matrix(EP[:vS_VRE_STOR]),
+        output_settings_d["WriteNSE"] ? Array{Float64, 3}(Array(value.(EP[:vNSE]))) : nothing,
+        output_settings_d["WriteEmissions"] ? _extract_output_matrix(EP[:eEmissionsByZone]) : nothing)
+end
+
+@doc raw"""
+    resource_time_scratch!(cache::OutputCache) -> Matrix{Float64}
+
+Clear and return the reusable `G x T` scratch matrix stored in `cache`.
+"""
+function resource_time_scratch!(cache::OutputCache)::Matrix{Float64}
+    fill!(cache.resource_time_scratch, 0.0)
+    return cache.resource_time_scratch
+end
+
+@doc raw"""
+    scaled_resource_time_matrix!(cache::OutputCache, data::Matrix{Float64})
+
+Return `data` scaled by `cache.scale_factor` while avoiding unnecessary
+allocations. If scaling is not needed, returns `data`; otherwise, writes the
+scaled result into cache scratch memory and returns that scratch matrix.
+"""
+function scaled_resource_time_matrix!(cache::OutputCache, data::Matrix{Float64})
+    if cache.scale_factor == 1
+        return data
+    end
+    scratch = resource_time_scratch!(cache)
+    copyto!(scratch, data)
+    rmul!(scratch, cache.scale_factor)
+    return scratch
+end
+
+@doc raw"""
+    materialize_output_blocks(blocks, block_ids::Vector{Vector{Int}}, T::Int)
+
+Stack a list of `blocks` (each with `T` columns) into one dense matrix and the
+corresponding flattened resource-id vector.
+
+# Arguments
+- `blocks`: Row blocks to stack.
+- `block_ids::Vector{Vector{Int}}`: Resource ids for each block row.
+- `T::Int`: Number of time steps (columns).
+
+# Returns
+- `(data, ids)`: `data::Matrix{Float64}` and `ids::Vector{Int}`.
+"""
+function materialize_output_blocks(blocks, block_ids::Vector{Vector{Int}}, T::Int)
+    total_rows = sum(length, block_ids; init = 0)
+    data = Matrix{Float64}(undef, total_rows, T)
+    ids = Vector{Int}(undef, total_rows)
+
+    next_row = 1
+    for (block, ids_block) in zip(blocks, block_ids)
+        row_count = length(ids_block)
+        row_range = next_row:(next_row + row_count - 1)
+        data[row_range, :] .= block
+        ids[row_range] .= ids_block
+        next_row += row_count
+    end
+
+    return data, ids
+end
\ No newline at end of file
diff --git a/src/write_outputs/reserves/write_rsv.jl b/src/write_outputs/reserves/write_rsv.jl
index 19ed31aadc..f1ad8b5457 100644
--- a/src/write_outputs/reserves/write_rsv.jl
+++ b/src/write_outputs/reserves/write_rsv.jl
@@ -30,8 +30,8 @@ function write_rsv(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
         rename!(total, auxNew_Names)
         rename!(unmet, auxNew_Names)
         dfRsv = vcat(dfRsv, unmet, total)
-        CSV.write(joinpath(path, "reserves.csv"),
-            dftranspose(dfRsv, false),
+        write_transposed_csv(joinpath(path, "reserves.csv"),
+            dfRsv,
             writeheader = false)
     end
 end
diff --git a/src/write_outputs/transmission/write_transmission_flows.jl b/src/write_outputs/transmission/write_transmission_flows.jl
index 5290d71afe..baf9cc56aa 100644
--- a/src/write_outputs/transmission/write_transmission_flows.jl
+++ b/src/write_outputs/transmission/write_transmission_flows.jl
@@ -22,7 +22,7 @@ function write_transmission_flows(path::AbstractString,
         dfFlow = hcat(dfFlow, DataFrame(flow, :auto))
         auxNew_Names = [Symbol("Line"); [Symbol("t$t") for t in 1:T]]
         rename!(dfFlow, auxNew_Names)
-        CSV.write(filepath, dftranspose(dfFlow, false), writeheader = false)
+        write_transposed_csv(filepath, dfFlow, writeheader = false)
 
         if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
             write_full_time_series_reconstruction(path, setup, dfFlow, "flow")
diff --git a/src/write_outputs/transmission/write_transmission_losses.jl b/src/write_outputs/transmission/write_transmission_losses.jl
index a76bca1180..c9b2394db5 100644
--- a/src/write_outputs/transmission/write_transmission_losses.jl
+++ b/src/write_outputs/transmission/write_transmission_losses.jl
@@ -27,8 +27,8 @@ function write_transmission_losses(path::AbstractString,
             auxNew_Names)
         total[:, 3:(T + 2)] .= sum(tlosses, dims = 1)
         dfTLosses = vcat(dfTLosses, total)
-        CSV.write(joinpath(path, "tlosses.csv"),
-            dftranspose(dfTLosses, false),
+        write_transposed_csv(joinpath(path, "tlosses.csv"),
+            dfTLosses,
             writeheader = false)
 
         if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
diff --git a/src/write_outputs/write_angles.jl b/src/write_outputs/write_angles.jl
index b93870354f..6c89dc358e 100644
--- a/src/write_outputs/write_angles.jl
+++ b/src/write_outputs/write_angles.jl
@@ -15,8 +15,6 @@ function write_angles(path::AbstractString, inputs::Dict, setup::Dict, EP::Model
     rename!(dfAngles, auxNew_Names)
 
     ## Linear configuration final output
-    CSV.write(joinpath(path, "angles.csv"),
-        dftranspose(dfAngles, false),
-        writeheader = false)
+    write_transposed_csv(joinpath(path, "angles.csv"), dfAngles, writeheader = false)
     return nothing
 end
diff --git a/src/write_outputs/write_capacityfactor.jl b/src/write_outputs/write_capacityfactor.jl
index c6f5352bf8..8fb308683c 100644
--- a/src/write_outputs/write_capacityfactor.jl
+++ b/src/write_outputs/write_capacityfactor.jl
@@ -1,10 +1,17 @@
 @doc raw"""
-	write_capacityfactor(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_capacityfactor(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing the capacity factor of different resources. For co-located VRE-storage resources, this
     value is calculated if the site has either or both a solar PV or wind resource.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_capacityfactor(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_capacityfactor(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]
     G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
     T = inputs["T"]     # Number of time steps (hours)
@@ -22,9 +29,9 @@ function write_capacityfactor(path::AbstractString, inputs::Dict, setup::Dict, E
         AnnualSum = zeros(G),
         Capacity = zeros(G),
         CapacityFactor = zeros(G))
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-    df.AnnualSum .= value.(EP[:vP]) * weight * scale_factor
-    df.Capacity .= value.(EP[:eTotalCap]) * scale_factor
+    scale_factor = cache.scale_factor
+    df.AnnualSum .= cache.vP * weight * scale_factor
+    df.Capacity .= cache.eTotalCap * scale_factor
 
     # The .data only works on DenseAxisArray variables or expressions
     # In contrast vP and eTotalCap are whole vectors / matrices
diff --git a/src/write_outputs/write_charge.jl b/src/write_outputs/write_charge.jl
index 35f2331847..07e3e66b3b 100644
--- a/src/write_outputs/write_charge.jl
+++ b/src/write_outputs/write_charge.jl
@@ -1,9 +1,16 @@
 @doc raw"""
-	write_charge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_charge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing the charging energy values of the different storage technologies.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_charge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_charge(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]   # Resources (objects) 
     resources = inputs["RESOURCE_NAMES"]    # Resource names
     zones = zone_id.(gen)
@@ -18,39 +25,38 @@ function write_charge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model
     FUSION = ids_with(gen, :fusion)
 
     weight = inputs["omega"]
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-
-    charge = Matrix[]
+    charge = Matrix{Float64}[]
     charge_ids = Vector{Int}[]
     if !isempty(STOR_ALL)
-        push!(charge, value.(EP[:vCHARGE]))
+        push!(charge, cache.vCHARGE)
         push!(charge_ids, STOR_ALL)
     end
     if !isempty(FLEX)
-        push!(charge, value.(EP[:vCHARGE_FLEX]))
+        push!(charge, cache.vCHARGE_FLEX)
         push!(charge_ids, FLEX)
     end
     if (setup["HydrogenMinimumProduction"] > 0) & (!isempty(ELECTROLYZER))
-        push!(charge, value.(EP[:vUSE]))
+        push!(charge, cache.vUSE)
         push!(charge_ids, ELECTROLYZER)
     end
     if !isempty(VS_STOR)
-        push!(charge, value.(EP[:vCHARGE_VRE_STOR]))
+        push!(charge, cache.vCHARGE_VRE_STOR)
         push!(charge_ids, VS_STOR)
     end
     if !isempty(FUSION)
         _, mat = prepare_fusion_parasitic_power(EP, inputs)
-        push!(charge, mat)
+        push!(charge, Matrix{Float64}(mat))
         push!(charge_ids, FUSION)
     end
     if !isempty(ALLAM_CYCLE_LOX)
-        push!(charge, value.(EP[:vCHARGE_ALLAM]))
+        push!(charge, cache.vCHARGE_ALLAM)
         push!(charge_ids, ALLAM_CYCLE_LOX)
     end
-    charge = reduce(vcat, charge, init = zeros(0, T))
-    charge_ids = reduce(vcat, charge_ids, init = Int[])
+    charge, charge_ids = materialize_output_blocks(charge, charge_ids, T)
 
-    charge *= scale_factor
+    if cache.scale_factor != 1
+        charge .*= cache.scale_factor
+    end
 
     df = DataFrame(Resource = resources[charge_ids],
         Zone = zones[charge_ids])
diff --git a/src/write_outputs/write_charging_cost.jl b/src/write_outputs/write_charging_cost.jl
index cabc4db135..2b5a249060 100644
--- a/src/write_outputs/write_charging_cost.jl
+++ b/src/write_outputs/write_charging_cost.jl
@@ -1,4 +1,17 @@
-function write_charging_cost(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+@doc raw"""
+	write_charging_cost(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+	    cache::OutputCache = build_output_cache(EP, inputs, setup))
+
+Function for writing charging cost from storage, flexible demand, electrolyzers,
+and co-located VRE-storage resources.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
+"""
+function write_charging_cost(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]  # Resources (objects)
     resources = inputs["RESOURCE_NAMES"] # Resource names
 
@@ -7,7 +20,6 @@ function write_charging_cost(path::AbstractString, inputs::Dict, setup::Dict, EP
     zones = zone_id.(gen)
 
     G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-    T = inputs["T"]     # Number of time steps (hours)
     STOR_ALL = inputs["STOR_ALL"]
     FLEX = inputs["FLEX"]
     ELECTROLYZER = inputs["ELECTROLYZER"]
@@ -16,32 +28,33 @@ function write_charging_cost(path::AbstractString, inputs::Dict, setup::Dict, EP
     FUSION = ids_with(gen, :fusion)
 
     weight = inputs["omega"]
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-
-    price = locational_marginal_price(EP, inputs, setup)
+    price = isnothing(cache.price) ? locational_marginal_price(EP, inputs, setup) :
+            cache.price
 
-    chargecost = zeros(G, T)
+    chargecost = resource_time_scratch!(cache)
     if !isempty(STOR_ALL)
-        chargecost[STOR_ALL, :] .= (value.(EP[:vCHARGE][STOR_ALL, :]).data) .*
+        chargecost[STOR_ALL, :] .= cache.vCHARGE .*
                                    transpose(price)[zone_id.(gen.Storage), :]
     end
     if !isempty(FLEX)
-        chargecost[FLEX, :] .= value.(EP[:vP][FLEX, :]) .*
+        chargecost[FLEX, :] .= cache.vP[FLEX, :] .*
                                transpose(price)[zone_id.(gen.FlexDemand), :]
     end
     if !isempty(ELECTROLYZER)
-        chargecost[ELECTROLYZER, :] .= (value.(EP[:vUSE][ELECTROLYZER, :]).data) .*
+        chargecost[ELECTROLYZER, :] .= cache.vUSE .*
                                        transpose(price)[zone_id.(gen.Electrolyzer), :]
     end
     if !isempty(VS_STOR)
-        chargecost[VS_STOR, :] .= value.(EP[:vCHARGE_VRE_STOR][VS_STOR, :].data) .*
+        chargecost[VS_STOR, :] .= cache.vCHARGE_VRE_STOR .*
                                   transpose(price)[zone_id.(gen[VS_STOR]), :]
     end
     if !isempty(FUSION)
         _, mat = prepare_fusion_parasitic_power(EP, inputs)
-        chargecost[FUSION, :] = mat
+        chargecost[FUSION, :] .= mat
+    end
+    if cache.scale_factor != 1
+        rmul!(chargecost, cache.scale_factor)
     end
-    chargecost *= scale_factor
 
     dfChargingcost = DataFrame(Region = regions,
         Resource = resources,
diff --git a/src/write_outputs/write_curtailment.jl b/src/write_outputs/write_curtailment.jl
index 4d41289750..a0a363f113 100644
--- a/src/write_outputs/write_curtailment.jl
+++ b/src/write_outputs/write_curtailment.jl
@@ -1,10 +1,17 @@
 @doc raw"""
-	write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing the curtailment values of the different variable renewable resources (both standalone and 
 	co-located).
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_curtailment(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]  # Resources (objects)
     resources = inputs["RESOURCE_NAMES"] # Resource names
     zones = zone_id.(gen)
@@ -15,11 +22,9 @@ function write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::
     VRE_STOR = inputs["VRE_STOR"]
 
     weight = inputs["omega"]
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-
-    curtailment = zeros(G, T)
-    curtailment[VRE, :] = (value.(EP[:eTotalCap][VRE]) .* inputs["pP_Max"][VRE, :] .-
-                           value.(EP[:vP][VRE, :]))
+    curtailment = resource_time_scratch!(cache)
+    curtailment[VRE, :] .= cache.eTotalCap[VRE] .* inputs["pP_Max"][VRE, :] .-
+                           cache.vP[VRE, :]
 
     if !isempty(VRE_STOR)
         SOLAR = setdiff(inputs["VS_SOLAR"], inputs["VS_WIND"])
@@ -48,7 +53,9 @@ function write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::
         end
     end
 
-    curtailment *= scale_factor
+    if cache.scale_factor != 1
+        rmul!(curtailment, cache.scale_factor)
+    end
 
     df = DataFrame(Resource = resources,
         Zone = zones,
diff --git a/src/write_outputs/write_emissions.jl b/src/write_outputs/write_emissions.jl
index 2e0c011f68..86a3fc9882 100644
--- a/src/write_outputs/write_emissions.jl
+++ b/src/write_outputs/write_emissions.jl
@@ -1,10 +1,17 @@
 @doc raw"""
-	write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for reporting time-dependent CO$_2$ emissions by zone.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 
 """
-function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_emissions(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     T = inputs["T"]     # Number of time steps (hours)
     Z = inputs["Z"]     # Number of zones
 
@@ -39,11 +46,10 @@ function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Mo
             dfEmissions = DataFrame(Zone = 1:Z, AnnualSum = Array{Float64}(undef, Z))
         end
 
-        emissions_by_zone = value.(EP[:eEmissionsByZone])
-        for i in 1:Z
-            dfEmissions[i, :AnnualSum] = sum(inputs["omega"] .* emissions_by_zone[i, :]) *
-                                         scale_factor
-        end
+        emissions_by_zone = isnothing(cache.eEmissionsByZone) ?
+                            Matrix{Float64}(value.(EP[:eEmissionsByZone])) :
+                            cache.eEmissionsByZone
+        dfEmissions[!, :AnnualSum] .= emissions_by_zone * inputs["omega"] * scale_factor
 
         if setup["WriteOutputs"] == "annual"
             total = DataFrame(["Total" sum(dfEmissions.AnnualSum)], [:Zone; :AnnualSum])
@@ -88,14 +94,16 @@ function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Mo
             end
             rename!(total, auxNew_Names)
             dfEmissions = vcat(dfEmissions, total)
-            CSV.write(joinpath(path, "emissions.csv"),
-                dftranspose(dfEmissions, false),
+            write_transposed_csv(joinpath(path, "emissions.csv"),
+                dfEmissions,
                 writeheader = false)
         end
         ## Aaron - Combined elseif setup["Dual_MIP"]==1 block with the first block since they were identical. Why do we have this third case? What is different about it?
     else
         # CO2 emissions by zone
-        emissions_by_zone = value.(EP[:eEmissionsByZone])
+        emissions_by_zone = isnothing(cache.eEmissionsByZone) ?
+                            Matrix{Float64}(value.(EP[:eEmissionsByZone])) :
+                            cache.eEmissionsByZone
         dfEmissions = hcat(DataFrame(Zone = 1:Z),
             DataFrame(AnnualSum = Array{Float64}(undef, Z)))
         for i in 1:Z
@@ -121,8 +129,8 @@ function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Mo
             end
             rename!(total, auxNew_Names)
             dfEmissions = vcat(dfEmissions, total)
-            CSV.write(joinpath(path, "emissions.csv"),
-                dftranspose(dfEmissions, false),
+            write_transposed_csv(joinpath(path, "emissions.csv"),
+                dfEmissions,
                 writeheader = false)
 
             if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
diff --git a/src/write_outputs/write_energy_revenue.jl b/src/write_outputs/write_energy_revenue.jl
index 3e0834bd1e..01154d6060 100644
--- a/src/write_outputs/write_energy_revenue.jl
+++ b/src/write_outputs/write_energy_revenue.jl
@@ -1,16 +1,22 @@
 @doc raw"""
-	write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing energy revenue from the different generation technologies.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_energy_revenue(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]
     regions = region.(gen)
     clusters = cluster.(gen)
     zones = zone_id.(gen)
 
     G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-    T = inputs["T"]     # Number of time steps (hours)
     FLEX = inputs["FLEX"]
     NONFLEX = setdiff(collect(1:G), FLEX)
     dfEnergyRevenue = DataFrame(Region = regions,
@@ -18,16 +24,17 @@ function write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, E
         Zone = zones,
         Cluster = clusters,
         AnnualSum = Array{Float64}(undef, G))
-    energyrevenue = zeros(G, T)
-    price = locational_marginal_price(EP, inputs, setup)
-    energyrevenue[NONFLEX, :] = value.(EP[:vP][NONFLEX, :]) .*
-                                transpose(price)[zone_id.(gen[NONFLEX]), :]
+    energyrevenue = resource_time_scratch!(cache)
+    price = isnothing(cache.price) ? locational_marginal_price(EP, inputs, setup) :
+            cache.price
+    energyrevenue[NONFLEX, :] .= cache.vP[NONFLEX, :] .*
+                                 transpose(price)[zone_id.(gen[NONFLEX]), :]
     if !isempty(FLEX)
-        energyrevenue[FLEX, :] = value.(EP[:vCHARGE_FLEX][FLEX, :]).data .*
-                                 transpose(price)[zone_id.(gen[FLEX]), :]
+        energyrevenue[FLEX, :] .= cache.vCHARGE_FLEX .*
+                                  transpose(price)[zone_id.(gen[FLEX]), :]
     end
-    if setup["ParameterScale"] == 1
-        energyrevenue *= ModelScalingFactor
+    if cache.scale_factor != 1
+        rmul!(energyrevenue, cache.scale_factor)
     end
     dfEnergyRevenue.AnnualSum .= energyrevenue * inputs["omega"]
     write_simple_csv(joinpath(path, "EnergyRevenue.csv"), dfEnergyRevenue)
diff --git a/src/write_outputs/write_fuel_consumption.jl b/src/write_outputs/write_fuel_consumption.jl
index 8385e3e5ce..a48892c82a 100644
--- a/src/write_outputs/write_fuel_consumption.jl
+++ b/src/write_outputs/write_fuel_consumption.jl
@@ -83,8 +83,9 @@ function write_fuel_consumption_ts(path::AbstractString,
     end
     dfPlantFuel_TS = hcat(dfPlantFuel_TS,
         DataFrame(tempts, [Symbol("t$t") for t in 1:T]))
-    CSV.write(joinpath(path, "FuelConsumption_plant_MMBTU.csv"),
-        dftranspose(dfPlantFuel_TS, false), header = false)
+    write_transposed_csv(joinpath(path, "FuelConsumption_plant_MMBTU.csv"),
+        dfPlantFuel_TS,
+        header = false)
 
     if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
         write_full_time_series_reconstruction(
diff --git a/src/write_outputs/write_nse.jl b/src/write_outputs/write_nse.jl
index 9d1c73e835..2823f7c4d7 100644
--- a/src/write_outputs/write_nse.jl
+++ b/src/write_outputs/write_nse.jl
@@ -1,9 +1,16 @@
 @doc raw"""
-	write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for reporting non-served energy for every model zone, time step and cost-segment.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_nse(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     T = inputs["T"]     # Number of time steps (hours)
     Z = inputs["Z"]     # Number of zones
     SEG = inputs["SEG"] # Number of demand curtailment segments
@@ -11,10 +18,12 @@ function write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
     dfNse = DataFrame(Segment = repeat(1:SEG, outer = Z),
         Zone = repeat(1:Z, inner = SEG),
         AnnualSum = zeros(SEG * Z))
-    nse = zeros(SEG * Z, T)
     scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
+    nse = zeros(SEG * Z, T)
+    nse_values = isnothing(cache.vNSE) ? Array{Float64, 3}(Array(value.(EP[:vNSE]))) :
+                 cache.vNSE
     for z in 1:Z
-        nse[((z - 1) * SEG + 1):(z * SEG), :] = value.(EP[:vNSE])[:, :, z] * scale_factor
+        nse[((z - 1) * SEG + 1):(z * SEG), :] .= nse_values[:, :, z] * scale_factor
     end
     dfNse.AnnualSum .= nse * inputs["omega"]
 
@@ -36,7 +45,7 @@ function write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
         rename!(total, auxNew_Names)
         dfNse = vcat(dfNse, total)
 
-        CSV.write(joinpath(path, "nse.csv"), dftranspose(dfNse, false), writeheader = false)
+        write_transposed_csv(joinpath(path, "nse.csv"), dfNse, writeheader = false)
 
         if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
             write_full_time_series_reconstruction(path, setup, dfNse, "nse")
diff --git a/src/write_outputs/write_outputs.jl b/src/write_outputs/write_outputs.jl
index 8779c3a783..15ad56fb2f 100644
--- a/src/write_outputs/write_outputs.jl
+++ b/src/write_outputs/write_outputs.jl
@@ -61,6 +61,8 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
         optimize!(EP)
     end
 
+    cache = build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
     if output_settings_d["WriteCosts"]
         elapsed_time_costs = @elapsed write_costs(path, inputs, setup, EP)
         println("Time elapsed for writing costs is")
@@ -74,37 +76,45 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
     end
 
     if output_settings_d["WritePower"] || output_settings_d["WriteNetRevenue"]
-        elapsed_time_power = @elapsed dfPower = write_power(path, inputs, setup, EP)
+        elapsed_time_power = @elapsed dfPower = write_power(path, inputs, setup, EP, cache)
         println("Time elapsed for writing power is")
         println(elapsed_time_power)
     end
 
     if output_settings_d["WriteCharge"]
-        elapsed_time_charge = @elapsed write_charge(path, inputs, setup, EP)
+        elapsed_time_charge = @elapsed write_charge(path, inputs, setup, EP, cache)
         println("Time elapsed for writing charge is")
         println(elapsed_time_charge)
     end
 
     if output_settings_d["WriteCapacityFactor"]
-        elapsed_time_capacityfactor = @elapsed write_capacityfactor(path, inputs, setup, EP)
+        elapsed_time_capacityfactor = @elapsed write_capacityfactor(path,
+            inputs,
+            setup,
+            EP,
+            cache)
         println("Time elapsed for writing capacity factor is")
         println(elapsed_time_capacityfactor)
     end
 
     if output_settings_d["WriteStorage"]
-        elapsed_time_storage = @elapsed write_storage(path, inputs, setup, EP)
+        elapsed_time_storage = @elapsed write_storage(path, inputs, setup, EP, cache)
         println("Time elapsed for writing storage is")
         println(elapsed_time_storage)
     end
 
     if output_settings_d["WriteCurtailment"]
-        elapsed_time_curtailment = @elapsed write_curtailment(path, inputs, setup, EP)
+        elapsed_time_curtailment = @elapsed write_curtailment(path,
+            inputs,
+            setup,
+            EP,
+            cache)
         println("Time elapsed for writing curtailment is")
         println(elapsed_time_curtailment)
     end
 
     if output_settings_d["WriteNSE"]
-        elapsed_time_nse = @elapsed write_nse(path, inputs, setup, EP)
+        elapsed_time_nse = @elapsed write_nse(path, inputs, setup, EP, cache)
         println("Time elapsed for writing nse is")
         println(elapsed_time_nse)
     end
@@ -139,7 +149,7 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
     end
 
     if output_settings_d["WriteEmissions"]
-        elapsed_time_emissions = @elapsed write_emissions(path, inputs, setup, EP)
+        elapsed_time_emissions = @elapsed write_emissions(path, inputs, setup, EP, cache)
         println("Time elapsed for writing emissions is")
         println(elapsed_time_emissions)
     end
@@ -174,7 +184,11 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
         end
         if !isempty(inputs["STOR_ALL"]) || !isempty(VS_STOR)
             if output_settings_d["WriteStorageDual"]
-                elapsed_time_stordual = @elapsed write_storagedual(path, inputs, setup, EP)
+                elapsed_time_stordual = @elapsed write_storagedual(path,
+                    inputs,
+                    setup,
+                    EP,
+                    cache)
                 println("Time elapsed for writing storage duals is")
                 println(elapsed_time_stordual)
             end
@@ -278,7 +292,7 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
         dfRegSubRevenue = DataFrame()
         if has_duals(EP) == 1
             if output_settings_d["WritePrice"]
-                elapsed_time_price = @elapsed write_price(path, inputs, setup, EP)
+                elapsed_time_price = @elapsed write_price(path, inputs, setup, EP, cache)
                 println("Time elapsed for writing price is")
                 println(elapsed_time_price)
             end
@@ -289,7 +303,8 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
                     path,
                     inputs,
                     setup,
-                    EP)
+                    EP,
+                    cache)
                 println("Time elapsed for writing energy revenue is")
                 println(elapsed_time_energy_rev)
             end
@@ -300,7 +315,8 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
                     path,
                     inputs,
                     setup,
-                    EP)
+                    EP,
+                    cache)
                 println("Time elapsed for writing charging cost is")
                 println(elapsed_time_charging_cost)
             end
@@ -537,18 +553,15 @@ function write_fulltimeseries(fullpath::AbstractString,
         dataOut::Matrix{Float64},
         dfOut::DataFrame)
     T = size(dataOut, 2)
-    dfOut = hcat(dfOut, DataFrame(dataOut, :auto))
-    auxNew_Names = [Symbol("Resource");
-                    Symbol("Zone");
-                    Symbol("AnnualSum");
-                    [Symbol("t$t") for t in 1:T]]
-    rename!(dfOut, auxNew_Names)
-    total = DataFrame(
-        ["Total" 0 sum(dfOut[!, :AnnualSum], init = 0.0) fill(0.0, (1, T))], auxNew_Names)
-    total[!, 4:(T + 3)] .= sum(dataOut, dims = 1, init = 0.0)
-    dfOut = vcat(dfOut, total)
-
-    CSV.write(fullpath, dftranspose(dfOut, false), writeheader = false)
+    for t in 1:T
+        dfOut[!, Symbol("t$t")] = dataOut[:, t]
+    end
+
+    total_row = Any["Total", 0, sum(dfOut[!, :AnnualSum], init = 0.0)]
+    append!(total_row, vec(sum(dataOut, dims = 1, init = 0.0)))
+    push!(dfOut, total_row)
+
+    write_transposed_csv(fullpath, dfOut, writeheader = false)
     return dfOut
 end
 
@@ -645,7 +658,10 @@ function write_full_time_series_reconstruction(
         path::AbstractString, setup::Dict, DF::DataFrame, name::String)
     FullTimeSeriesFolder = setup["OutputFullTimeSeriesFolder"]
     output_path = joinpath(path, FullTimeSeriesFolder)
-    dfOut_full = full_time_series_reconstruction(path, setup, dftranspose(DF, false))
+    dfOut_full = full_time_series_reconstruction(
+        path,
+        setup,
+        transpose_output_dataframe(DF))
     CSV.write(joinpath(output_path, "$name.csv"), dfOut_full, header = false)
     return nothing
 end
diff --git a/src/write_outputs/write_power.jl b/src/write_outputs/write_power.jl
index 995a2d941b..0f69d23792 100644
--- a/src/write_outputs/write_power.jl
+++ b/src/write_outputs/write_power.jl
@@ -1,9 +1,16 @@
 @doc raw"""
-	write_power(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_power(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing the different values of power generated by the different technologies in operation.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_power(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_power(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]   # Resources (objects)
     resources = inputs["RESOURCE_NAMES"]    # Resource names
     zones = zone_id.(gen)
@@ -12,11 +19,7 @@ function write_power(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
     T = inputs["T"]     # Number of time steps (hours)
     
     weight = inputs["omega"]
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-
-    # Power injected by each resource in each time step
-    power = value.(EP[:vP])
-    power *= scale_factor
+    power = scaled_resource_time_matrix!(cache, cache.vP)
 
     df = DataFrame(Resource = resources,
         Zone = zones,
diff --git a/src/write_outputs/write_power_balance.jl b/src/write_outputs/write_power_balance.jl
index 69dd596bc5..55b2e51044 100644
--- a/src/write_outputs/write_power_balance.jl
+++ b/src/write_outputs/write_power_balance.jl
@@ -117,8 +117,8 @@ function write_power_balance(path::AbstractString, inputs::Dict, setup::Dict, EP
                         Symbol("AnnualSum");
                         [Symbol("t$t") for t in 1:T]]
         rename!(dfPowerBalance, auxNew_Names)
-        CSV.write(joinpath(path, "power_balance.csv"),
-            dftranspose(dfPowerBalance, false),
+        write_transposed_csv(joinpath(path, "power_balance.csv"),
+            dfPowerBalance,
             writeheader = false)
 
         if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
diff --git a/src/write_outputs/write_price.jl b/src/write_outputs/write_price.jl
index bda7c05a16..0f8c5f3328 100644
--- a/src/write_outputs/write_price.jl
+++ b/src/write_outputs/write_price.jl
@@ -1,9 +1,16 @@
 @doc raw"""
-	write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for reporting marginal electricity price for each model zone and time step. Marginal electricity price is equal to the dual variable of the power balance constraint. If GenX is configured as a mixed integer linear program, then this output is only generated if `WriteShadowPrices` flag is activated. If configured as a linear program (i.e. linearized unit commitment or economic dispatch) then output automatically available.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_price(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     T = inputs["T"]     # Number of time steps (hours)
     Z = inputs["Z"]     # Number of zones
 
@@ -11,16 +18,15 @@ function write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
     # Electricity price: Dual variable of hourly power balance constraint = hourly price
     dfPrice = DataFrame(Zone = 1:Z) # The unit is $/MWh
     # Dividing dual variable for each hour with corresponding hourly weight to retrieve marginal cost of generation
-    price = locational_marginal_price(EP, inputs, setup)
+        price = isnothing(cache.price) ? locational_marginal_price(EP, inputs, setup) :
+            cache.price
     dfPrice = hcat(dfPrice, DataFrame(transpose(price), :auto))
 
     auxNew_Names = [Symbol("Zone"); [Symbol("t$t") for t in 1:T]]
     rename!(dfPrice, auxNew_Names)
 
     ## Linear configuration final output
-    CSV.write(joinpath(path, "prices.csv"),
-        dftranspose(dfPrice, false),
-        writeheader = false)
+    write_transposed_csv(joinpath(path, "prices.csv"), dfPrice, writeheader = false)
 
     if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
         write_full_time_series_reconstruction(path, setup, dfPrice, "prices")
diff --git a/src/write_outputs/write_reliability.jl b/src/write_outputs/write_reliability.jl
index afb3a3284c..0a42fcfc71 100644
--- a/src/write_outputs/write_reliability.jl
+++ b/src/write_outputs/write_reliability.jl
@@ -17,8 +17,8 @@ function write_reliability(path::AbstractString, inputs::Dict, setup::Dict, EP::
     auxNew_Names = [Symbol("Zone"); [Symbol("t$t") for t in 1:T]]
     rename!(dfReliability, auxNew_Names)
 
-    CSV.write(joinpath(path, "reliability.csv"),
-        dftranspose(dfReliability, false),
+    write_transposed_csv(joinpath(path, "reliability.csv"),
+        dfReliability,
         header = false)
 
     if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
diff --git a/src/write_outputs/write_storage.jl b/src/write_outputs/write_storage.jl
index 4e64c0f9bd..c3af09ed84 100644
--- a/src/write_outputs/write_storage.jl
+++ b/src/write_outputs/write_storage.jl
@@ -1,15 +1,21 @@
 @doc raw"""
-	write_storage(path::AbstractString, inputs::Dict,setup::Dict, EP::Model)
+    write_storage(path::AbstractString, inputs::Dict,setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for writing the capacities of different storage technologies, including hydro reservoir, flexible storage tech etc.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_storage(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_storage(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]   # Resources (objects)
     resources = inputs["RESOURCE_NAMES"]   # Resource names
     zones = zone_id.(gen)
 
     T = inputs["T"]     # Number of time steps (hours)
-    G = inputs["G"]
     STOR_ALL = inputs["STOR_ALL"]
     HYDRO_RES = inputs["HYDRO_RES"]
     FLEX = inputs["FLEX"]
@@ -17,25 +23,29 @@ function write_storage(path::AbstractString, inputs::Dict, setup::Dict, EP::Mode
     VS_STOR = !isempty(VRE_STOR) ? inputs["VS_STOR"] : []
 
     weight = inputs["omega"]
-    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-
-    stored = Matrix[]
+    stored = Matrix{Float64}[]
+    stored_groups = Vector{Int}[]
     if !isempty(STOR_ALL)
-        push!(stored, value.(EP[:vS]))
+        push!(stored, cache.vS)
+        push!(stored_groups, STOR_ALL)
     end
     if !isempty(HYDRO_RES)
-        push!(stored, value.(EP[:vS_HYDRO]))
+        push!(stored, cache.vS_HYDRO)
+        push!(stored_groups, HYDRO_RES)
     end
     if !isempty(FLEX)
-        push!(stored, value.(EP[:vS_FLEX]))
+        push!(stored, cache.vS_FLEX)
+        push!(stored_groups, FLEX)
     end
     if !isempty(VS_STOR)
-        push!(stored, value.(EP[:vS_VRE_STOR]))
+        push!(stored, cache.vS_VRE_STOR)
+        push!(stored_groups, VS_STOR)
+    end
+    stored, stored_ids = materialize_output_blocks(stored, stored_groups, T)
+    if cache.scale_factor != 1
+        stored .*= cache.scale_factor
     end
-    stored = reduce(vcat, stored, init = zeros(0, T))
-    stored *= scale_factor
 
-    stored_ids = convert(Vector{Int}, vcat(STOR_ALL, HYDRO_RES, FLEX, VS_STOR))
     df = DataFrame(Resource = resources[stored_ids],
         Zone = zones[stored_ids])
     df.AnnualSum = stored * weight
diff --git a/src/write_outputs/write_storagedual.jl b/src/write_outputs/write_storagedual.jl
index 90eae94ff2..1b5638b923 100644
--- a/src/write_outputs/write_storagedual.jl
+++ b/src/write_outputs/write_storagedual.jl
@@ -1,9 +1,16 @@
 @doc raw"""
-	write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+    write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model,
+        cache::OutputCache = build_output_cache(EP, inputs, setup))
 
 Function for reporting dual of storage level (state of charge) balance of each resource in each time step.
+The optional `cache` argument allows callers to reuse extracted model outputs across
+multiple write functions to reduce memory allocations.
 """
-function write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+function write_storagedual(path::AbstractString,
+    inputs::Dict,
+    setup::Dict,
+    EP::Model,
+    cache::OutputCache = build_output_cache(EP, inputs, setup))
     gen = inputs["RESOURCES"]
     zones = zone_id.(gen)
 
@@ -23,25 +30,25 @@ function write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::
 
     # # Dual of storage level (state of charge) balance of each resource in each time step
     dfStorageDual = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
-    dual_values = zeros(G, T)
+    dual_values = resource_time_scratch!(cache)
 
     # Loop over W separately hours_per_subperiod
     if !isempty(STOR_ALL)
         STOR_ALL_NONLDS = setdiff(STOR_ALL, inputs["STOR_LONG_DURATION"])
         STOR_ALL_LDS = intersect(STOR_ALL, inputs["STOR_LONG_DURATION"])
-        dual_values[STOR_ALL, INTERIOR_SUBPERIODS] = (dual.(EP[:cSoCBalInterior][
+        dual_values[STOR_ALL, INTERIOR_SUBPERIODS] .= (dual.(EP[:cSoCBalInterior][
             INTERIOR_SUBPERIODS,
             STOR_ALL]).data ./ inputs["omega"][INTERIOR_SUBPERIODS])'
-        dual_values[STOR_ALL_NONLDS, START_SUBPERIODS] = (dual.(EP[:cSoCBalStart][
+        dual_values[STOR_ALL_NONLDS, START_SUBPERIODS] .= (dual.(EP[:cSoCBalStart][
             START_SUBPERIODS,
             STOR_ALL_NONLDS]).data ./ inputs["omega"][START_SUBPERIODS])'
         if !isempty(STOR_ALL_LDS)
             if inputs["REP_PERIOD"] > 1
-                dual_values[STOR_ALL_LDS, START_SUBPERIODS] = (dual.(EP[:cSoCBalLongDurationStorageStart][
+                dual_values[STOR_ALL_LDS, START_SUBPERIODS] .= (dual.(EP[:cSoCBalLongDurationStorageStart][
                     1:REP_PERIOD,
                     STOR_ALL_LDS]).data ./ inputs["omega"][START_SUBPERIODS])'
             else
-                dual_values[STOR_ALL_LDS, START_SUBPERIODS] = (dual.(EP[:cSoCBalStart][
+                dual_values[STOR_ALL_LDS, START_SUBPERIODS] .= (dual.(EP[:cSoCBalStart][
                     START_SUBPERIODS,
                     STOR_ALL_LDS]).data ./ inputs["omega"][START_SUBPERIODS])'
             end
@@ -49,19 +56,19 @@ function write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::
     end
 
     if !isempty(VRE_STOR)
-        dual_values[VS_STOR, INTERIOR_SUBPERIODS] = ((dual.(EP[:cSoCBalInterior_VRE_STOR][
+        dual_values[VS_STOR, INTERIOR_SUBPERIODS] .= ((dual.(EP[:cSoCBalInterior_VRE_STOR][
             VS_STOR,
             INTERIOR_SUBPERIODS]).data)' ./ inputs["omega"][INTERIOR_SUBPERIODS])'
-        dual_values[VS_NONLDS, START_SUBPERIODS] = ((dual.(EP[:cSoCBalStart_VRE_STOR][
+        dual_values[VS_NONLDS, START_SUBPERIODS] .= ((dual.(EP[:cSoCBalStart_VRE_STOR][
             VS_NONLDS,
             START_SUBPERIODS]).data)' ./ inputs["omega"][START_SUBPERIODS])'
         if !isempty(VS_LDS)
             if inputs["REP_PERIOD"] > 1
-                dual_values[VS_LDS, START_SUBPERIODS] = ((dual.(EP[:cVreStorSoCBalLongDurationStorageStart][
+                dual_values[VS_LDS, START_SUBPERIODS] .= ((dual.(EP[:cVreStorSoCBalLongDurationStorageStart][
                     VS_LDS,
                     1:REP_PERIOD]).data)' ./ inputs["omega"][START_SUBPERIODS])'
             else
-                dual_values[VS_LDS, START_SUBPERIODS] = ((dual.(EP[:cSoCBalStart_VRE_STOR][
+                dual_values[VS_LDS, START_SUBPERIODS] .= ((dual.(EP[:cSoCBalStart_VRE_STOR][
                     VS_LDS,
                     START_SUBPERIODS]).data)' ./ inputs["omega"][START_SUBPERIODS])'
             end
@@ -76,8 +83,8 @@ function write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::
     rename!(dfStorageDual,
         [Symbol("Resource"); Symbol("Zone"); [Symbol("t$t") for t in 1:T]])
 
-    CSV.write(joinpath(path, "storagebal_duals.csv"),
-        dftranspose(dfStorageDual, false),
+    write_transposed_csv(joinpath(path, "storagebal_duals.csv"),
+        dfStorageDual,
         header = false)
 
     if setup["OutputFullTimeSeries"] == 1 && setup["TimeDomainReduction"] == 1
diff --git a/test/runtests.jl b/test/runtests.jl
index 57b4d68454..f9e7216e39 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -59,6 +59,10 @@ end
     @testset "Retrofit" begin
         include("test_retrofit.jl")
     end
+
+    @testset "Output cache and transpose" begin
+        include("test_output_cache.jl")
+    end
 end
 
 # Test writing outputs
diff --git a/test/test_output_cache.jl b/test/test_output_cache.jl
new file mode 100644
index 0000000000..24f5d3bc2d
--- /dev/null
+++ b/test/test_output_cache.jl
@@ -0,0 +1,260 @@
+using GenX, JuMP, HiGHS, CSV, DataFrames, Test
+
+function make_test_model(G, T, Z)
+    opt = optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false)
+    EP = Model(opt)
+
+    @variable(EP, vP[1:G, 1:T] >= 0)
+    @variable(EP, eTotalCap[1:G] >= 0)
+    @variable(EP, vNSE[1:1, 1:T, 1:Z] >= 0)
+    @variable(EP, eEmissionsByZone[1:Z, 1:T] >= 0)
+
+    @constraint(EP, [i = 1:G, t = 1:T], vP[i, t] == Float64(i * t))
+    @constraint(EP, [i = 1:G], eTotalCap[i] == Float64(i))
+    @constraint(EP, [s = 1:1, t = 1:T, z = 1:Z], vNSE[s, t, z] == 0.0)
+    @constraint(EP, [z = 1:Z, t = 1:T], eEmissionsByZone[z, t] == 0.0)
+
+    optimize!(EP)
+    return EP
+end
+
+function make_test_inputs(G, T)
+    return Dict(
+        "G" => G,
+        "T" => T,
+        "STOR_ALL" => Int[],
+        "FLEX" => Int[],
+        "ELECTROLYZER" => Int[],
+        "VRE_STOR" => Int[],
+        "ALLAM_CYCLE_LOX" => Int[],
+        "HYDRO_RES" => Int[],
+    )
+end
+
+function make_output_settings(overrides...)
+    d = Dict{String, Any}(
+        "WritePower" => false,
+        "WriteCapacityFactor" => false,
+        "WriteCurtailment" => false,
+        "WriteCharge" => false,
+        "WriteChargingCost" => false,
+        "WriteNetRevenue" => false,
+        "WriteStorage" => false,
+        "WriteStorageDual" => false,
+        "WriteEnergyRevenue" => false,
+        "WritePrice" => false,
+        "WriteNSE" => false,
+        "WriteEmissions" => false,
+    )
+    for (key, value) in overrides
+        d[key] = value
+    end
+    return d
+end
+
+function make_setup(output_settings_d)
+    return Dict(
+        "ParameterScale" => 0,
+        "WriteShadowPrices" => 0,
+        "WriteOutputsSettingsDict" => output_settings_d,
+    )
+end
+
+@testset "Output cache and transpose" begin
+    @testset "transpose_output_dataframe" begin
+        @testset "basic correctness and shape" begin
+            df = DataFrame(A = [1, 2, 3], B = [4, 5, 6], C = [7, 8, 9])
+            result = GenX.transpose_output_dataframe(df)
+
+            @test nrow(result) == ncol(df)
+            @test ncol(result) == nrow(df) + 1
+            @test result[!, :Row] == names(df)
+            @test result[1, :x1] == 1
+            @test result[2, :x1] == 4
+            @test result[3, :x1] == 7
+            @test result[1, :x2] == 2
+            @test result[1, :x3] == 3
+        end
+
+        @testset "withhead=true uses first column values as column names" begin
+            df = DataFrame(label = ["alpha", "beta", "gamma"], v1 = [10, 20, 30], v2 = [100, 200, 300])
+            result = GenX.transpose_output_dataframe(df; withhead = true)
+
+            @test collect(propertynames(result)) == [:Row, :alpha, :beta, :gamma]
+            @test result[1, :alpha] == "alpha"
+            @test result[2, :alpha] == 10
+            @test result[3, :alpha] == 100
+        end
+
+        @testset "input DataFrame is not mutated" begin
+            df = DataFrame(X = [1, 2], Y = [3, 4])
+            original_nrow = nrow(df)
+            original_ncol = ncol(df)
+            original_values = copy(Matrix(df))
+
+            GenX.transpose_output_dataframe(df)
+
+            @test nrow(df) == original_nrow
+            @test ncol(df) == original_ncol
+            @test Matrix(df) == original_values
+        end
+
+        @testset "single-row DataFrame transposes correctly" begin
+            df = DataFrame(p = [42], q = [99], r = [7])
+            result = GenX.transpose_output_dataframe(df)
+
+            @test nrow(result) == 3
+            @test ncol(result) == 2
+            @test result[!, :Row] == ["p", "q", "r"]
+            @test result[1, :x1] == 42
+            @test result[2, :x1] == 99
+            @test result[3, :x1] == 7
+        end
+    end
+
+    @testset "write_transposed_csv" begin
+        @testset "round-trip file shape and first column" begin
+            df = DataFrame(A = [1, 2, 3], B = [4, 5, 6])
+
+            mktempdir() do dir
+                path = joinpath(dir, "out.csv")
+                GenX.write_transposed_csv(path, df, writeheader = false)
+
+                on_disk = CSV.read(path, DataFrame; header = false)
+
+                @test nrow(on_disk) == ncol(df)
+                @test ncol(on_disk) == nrow(df) + 1
+                @test on_disk[!, 1] == string.(names(df))
+            end
+        end
+    end
+
+    @testset "build_output_cache" begin
+        G, T, Z = 4, 3, 2
+        EP = make_test_model(G, T, Z)
+
+        @testset "non-selective (backward compat)" begin
+            output_settings_d = make_output_settings()
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = false)
+
+            @test size(cache.vP) == (G, T)
+            @test cache.vP[2, 3] ≈ 6.0
+            @test length(cache.eTotalCap) == G
+            @test cache.eTotalCap[3] ≈ 3.0
+            @test isnothing(cache.vNSE)
+            @test isnothing(cache.eEmissionsByZone)
+        end
+
+        @testset "selective - all flags false" begin
+            output_settings_d = make_output_settings()
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
+            @test size(cache.vP) == (0, 0)
+            @test length(cache.eTotalCap) == 0
+            @test size(cache.resource_time_scratch) == (0, 0)
+            @test isnothing(cache.vNSE)
+            @test isnothing(cache.eEmissionsByZone)
+            @test isnothing(cache.vCHARGE)
+            @test isnothing(cache.vS)
+        end
+
+        @testset "selective - only WritePower=true" begin
+            output_settings_d = make_output_settings("WritePower" => true)
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
+            @test size(cache.vP) == (G, T)
+            @test cache.vP[1, 1] ≈ 1.0
+            @test length(cache.eTotalCap) == 0
+            @test isnothing(cache.vNSE)
+            @test isnothing(cache.eEmissionsByZone)
+        end
+
+        @testset "selective - only WriteNSE=true" begin
+            output_settings_d = make_output_settings("WriteNSE" => true)
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
+            @test !isnothing(cache.vNSE)
+            @test size(cache.vNSE) == (1, T, Z)
+            @test all(==(0.0), cache.vNSE)
+            @test size(cache.vP) == (0, 0)
+        end
+
+        @testset "selective - only WriteEmissions=true" begin
+            output_settings_d = make_output_settings("WriteEmissions" => true)
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
+            @test !isnothing(cache.eEmissionsByZone)
+            @test size(cache.eEmissionsByZone) == (Z, T)
+            @test all(==(0.0), cache.eEmissionsByZone)
+            @test size(cache.vP) == (0, 0)
+        end
+
+        @testset "selective - WriteCapacityFactor=true pulls both vP and eTotalCap" begin
+            output_settings_d = make_output_settings("WriteCapacityFactor" => true)
+            setup = make_setup(output_settings_d)
+            inputs = make_test_inputs(G, T)
+
+            cache = GenX.build_output_cache(EP, inputs, setup, output_settings_d; selective = true)
+
+            @test size(cache.vP) == (G, T)
+            @test length(cache.eTotalCap) == G
+            @test cache.eTotalCap[4] ≈ 4.0
+        end
+    end
+
+    @testset "resource_time_scratch! and scaled_resource_time_matrix!" begin
+        function make_cache(scale_factor, scratch_val = 0.0)
+            scratch = fill(scratch_val, 4, 3)
+            return GenX.OutputCache(
+                scale_factor,
+                scratch,
+                nothing,
+                zeros(Float64, 4, 3),
+                zeros(Float64, 4),
+                nothing, nothing, nothing,
+                nothing, nothing,
+                nothing, nothing, nothing,
+                nothing,
+                nothing,
+                nothing,
+            )
+        end
+
+        @testset "resource_time_scratch! zeroes the scratch matrix" begin
+            cache = make_cache(1.0, 99.0)
+            result = GenX.resource_time_scratch!(cache)
+            @test all(==(0.0), result)
+            @test result === cache.resource_time_scratch
+        end
+
+        @testset "scaled_resource_time_matrix! with scale_factor == 1 returns identity" begin
+            cache = make_cache(1.0)
+            data = rand(Float64, 4, 3)
+            result = GenX.scaled_resource_time_matrix!(cache, data)
+            @test result === data
+        end
+
+        @testset "scaled_resource_time_matrix! with scale_factor == 2.0 scales correctly" begin
+            cache = make_cache(2.0)
+            data = ones(Float64, 4, 3)
+            result = GenX.scaled_resource_time_matrix!(cache, data)
+            @test all(==(2.0), result)
+            @test result !== data
+        end
+    end
+end