diff --git a/docs/src/format.md b/docs/src/format.md
index 20b0a99..b43194a 100644
--- a/docs/src/format.md
+++ b/docs/src/format.md
@@ -9,6 +9,7 @@ An instance of the stochastic security-constrained unit commitment (SCUC) proble
 * [Parameters](#Parameters)
 * [Buses](#Buses)
 * [Generators](#Generators)
+* [Storage units](#Storage-units)
 * [Price-sensitive loads](#Price-sensitive-loads)
 * [Transmission lines](#Transmission-lines)
 * [Reserves](#Reserves)
@@ -173,6 +174,81 @@ Note that this curve also specifies the production limits. Specifically, the fir
 }
 ```
 
+### Storage units
+
+This section describes energy storage units in the system which charge and discharge power. The storage units consume power while charging, and generate power while discharging. 
+
+| Key               | Description                                       | Default  | Time series? | Uncertain?
+| :---------------- | :------------------------------------------------ | :------: | :------------: | :----: 
+| `Bus`             | Bus where the storage unit is located. Multiple storage units may be placed at the same bus. | Required | No | Yes
+| `Minimum level (MWh)`  | Minimum of energy level this storage unit may contain. | `0.0` | Yes | Yes
+| `Maximum level (MWh)`  | Maximum of energy level this storage unit may contain. | Required | Yes | Yes
+| `Allow simultaneous charging and discharging`  | If `false`, the storage unit is not allowed to charge and discharge at the same time (Boolean). | `true` | Yes | Yes
+| `Charge cost ($/MW)`  | Cost incurred for charging each MW of power into this storage unit. | Required | Yes | Yes
+| `Discharge cost ($/MW)`  | Cost incurred for discharging each MW of power from this storage unit. | Required | Yes | Yes
+| `Charge efficiency`  | Efficiency rate to charge power into this storage unit. This value must be greater than or equal to `0.0`, and less than or equal to `1.0`. | `1.0` | Yes | Yes
+| `Discharge efficiency`  | Efficiency rate to discharge power from this storage unit. This value must be greater than or equal to `0.0`, and less than or equal to `1.0`. | `1.0` | Yes | Yes
+| `Loss factor`  | The energy dissipation rate of this storage unit. This value must be greater than or equal to `0.0`, and less than or equal to `1.0`. | `0.0` | Yes | Yes
+| `Minimum charge rate (MW)`  | Minimum amount of power rate this storage unit may charge. | `0.0` | Yes | Yes
+| `Maximum charge rate (MW)`  | Maximum amount of power rate this storage unit may charge. | Required | Yes | Yes
+| `Minimum discharge rate (MW)`  | Minimum amount of power rate this storage unit may discharge. | `0.0` | Yes | Yes
+| `Maximum discharge rate (MW)`  | Maximum amount of power rate this storage unit may discharge. | Required | Yes | Yes
+| `Initial level (MWh)`  | Amount of energy this storage unit at time step `-1`, immediately before the planning horizon starts. | `0.0` | No | Yes
+| `Last period minimum level (MWh)`  | Minimum of energy level this storage unit may contain in the last time step. By default, this value is the same as the last value of `Minimum level (MWh)`. | `Minimum level (MWh)` | No | Yes
+| `Last period maximum level (MWh)`  | Maximum of energy level this storage unit may contain in the last time step. By default, this value is the same as the last value of `Maximum level (MWh)`. | `Maximum level (MWh)` | No | Yes
+
+#### Example
+```json
+{
+    "Storage units": {
+        "su1": {
+            "Bus": "b2",
+            "Maximum level (MWh)": 100.0,
+            "Charge cost ($/MW)": 2.0,
+            "Discharge cost ($/MW)": 2.5,
+            "Maximum charge rate (MW)": 10.0,
+            "Maximum discharge rate (MW)": 8.0
+        },
+        "su2": {
+            "Bus": "b2",
+            "Minimum level (MWh)": 10.0,
+            "Maximum level (MWh)": 100.0,
+            "Allow simultaneous charging and discharging": false,
+            "Charge cost ($/MW)": 3.0,
+            "Discharge cost ($/MW)": 3.5,
+            "Charge efficiency": 0.8,
+            "Discharge efficiency": 0.85,
+            "Loss factor": 0.01,
+            "Minimum charge rate (MW)": 5.0,
+            "Maximum charge rate (MW)": 10.0,
+            "Minimum discharge rate (MW)": 2.0,
+            "Maximum discharge rate (MW)": 10.0,
+            "Initial level (MWh)": 70.0,
+            "Last period minimum level (MWh)": 80.0,
+            "Last period maximum level (MWh)": 85.0
+        },
+        "su3": {
+            "Bus": "b9",
+            "Minimum level (MWh)": [10.0, 11.0, 12.0, 13.0],
+            "Maximum level (MWh)": [100.0, 110.0, 120.0, 130.0],
+            "Allow simultaneous charging and discharging": [false, false, true, true],
+            "Charge cost ($/MW)": [2.0, 2.1, 2.2, 2.3],
+            "Discharge cost ($/MW)": [1.0, 1.1, 1.2, 1.3],
+            "Charge efficiency": [0.8, 0.81, 0.82, 0.82],
+            "Discharge efficiency": [0.85, 0.86, 0.87, 0.88],
+            "Loss factor": [0.01, 0.01, 0.02, 0.02],
+            "Minimum charge rate (MW)": [5.0, 5.1, 5.2, 5.3],
+            "Maximum charge rate (MW)": [10.0, 10.1, 10.2, 10.3],
+            "Minimum discharge rate (MW)": [4.0, 4.1, 4.2, 4.3],
+            "Maximum discharge rate (MW)": [8.0, 8.1, 8.2, 8.3],
+            "Initial level (MWh)": 20.0,
+            "Last period minimum level (MWh)": 21.0,
+            "Last period maximum level (MWh)": 22.0
+        }
+    }
+}
+```
+
 ### Price-sensitive loads
 
 This section describes components in the system which may increase or reduce their energy consumption according to the energy prices. Fixed loads (as described in the `buses` section) are always served, regardless of the price, unless there is significant congestion in the system or insufficient production capacity. Price-sensitive loads, on the other hand, are only served if it is economical to do so. 
diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl
index 033c6df..d46d71e 100644
--- a/src/UnitCommitment.jl
+++ b/src/UnitCommitment.jl
@@ -36,6 +36,7 @@ include("model/formulations/base/sensitivity.jl")
 include("model/formulations/base/system.jl")
 include("model/formulations/base/unit.jl")
 include("model/formulations/base/punit.jl")
+include("model/formulations/base/storage.jl")
 include("model/formulations/CarArr2006/pwlcosts.jl")
 include("model/formulations/DamKucRajAta2016/ramp.jl")
 include("model/formulations/Gar1962/pwlcosts.jl")
diff --git a/src/instance/read.jl b/src/instance/read.jl
index 8c33657..10365ab 100644
--- a/src/instance/read.jl
+++ b/src/instance/read.jl
@@ -136,6 +136,7 @@ function _from_json(json; repair = true)::UnitCommitmentScenario
     loads = PriceSensitiveLoad[]
     reserves = Reserve[]
     profiled_units = ProfiledUnit[]
+    storage_units = StorageUnit[]
 
     function scalar(x; default = nothing)
         x !== nothing || return default
@@ -196,6 +197,7 @@ function _from_json(json; repair = true)::UnitCommitmentScenario
             ThermalUnit[],
             PriceSensitiveLoad[],
             ProfiledUnit[],
+            StorageUnit[],
         )
         name_to_bus[bus_name] = bus
         push!(buses, bus)
@@ -405,6 +407,52 @@ function _from_json(json; repair = true)::UnitCommitmentScenario
         end
     end
 
+    # Read storage units 
+    if "Storage units" in keys(json)
+        for (storage_name, dict) in json["Storage units"]
+            bus = name_to_bus[dict["Bus"]]
+            min_level =
+                timeseries(scalar(dict["Minimum level (MWh)"], default = 0.0))
+            max_level = timeseries(dict["Maximum level (MWh)"])
+            storage = StorageUnit(
+                storage_name,
+                bus,
+                min_level,
+                max_level,
+                timeseries(
+                    scalar(
+                        dict["Allow simultaneous charging and discharging"],
+                        default = true,
+                    ),
+                ),
+                timeseries(dict["Charge cost (\$/MW)"]),
+                timeseries(dict["Discharge cost (\$/MW)"]),
+                timeseries(scalar(dict["Charge efficiency"], default = 1.0)),
+                timeseries(scalar(dict["Discharge efficiency"], default = 1.0)),
+                timeseries(scalar(dict["Loss factor"], default = 0.0)),
+                timeseries(
+                    scalar(dict["Minimum charge rate (MW)"], default = 0.0),
+                ),
+                timeseries(dict["Maximum charge rate (MW)"]),
+                timeseries(
+                    scalar(dict["Minimum discharge rate (MW)"], default = 0.0),
+                ),
+                timeseries(dict["Maximum discharge rate (MW)"]),
+                scalar(dict["Initial level (MWh)"], default = 0.0),
+                scalar(
+                    dict["Last period minimum level (MWh)"],
+                    default = min_level[T],
+                ),
+                scalar(
+                    dict["Last period maximum level (MWh)"],
+                    default = max_level[T],
+                ),
+            )
+            push!(bus.storage_units, storage)
+            push!(storage_units, storage)
+        end
+    end
+
     scenario = UnitCommitmentScenario(
         name = scenario_name,
         probability = probability,
@@ -425,6 +473,8 @@ function _from_json(json; repair = true)::UnitCommitmentScenario
         thermal_units = thermal_units,
         profiled_units_by_name = Dict(pu.name => pu for pu in profiled_units),
         profiled_units = profiled_units,
+        storage_units_by_name = Dict(su.name => su for su in storage_units),
+        storage_units = storage_units,
         isf = spzeros(Float64, length(lines), length(buses) - 1),
         lodf = spzeros(Float64, length(lines), length(lines)),
     )
diff --git a/src/instance/structs.jl b/src/instance/structs.jl
index 39725af..3c1b635 100644
--- a/src/instance/structs.jl
+++ b/src/instance/structs.jl
@@ -9,6 +9,7 @@ mutable struct Bus
     thermal_units::Vector
     price_sensitive_loads::Vector
     profiled_units::Vector
+    storage_units::Vector
 end
 
 mutable struct CostSegment
@@ -83,6 +84,26 @@ mutable struct ProfiledUnit
     cost::Vector{Float64}
 end
 
+mutable struct StorageUnit
+    name::String
+    bus::Bus
+    min_level::Vector{Float64}
+    max_level::Vector{Float64}
+    simultaneous_charge_and_discharge::Vector{Bool}
+    charge_cost::Vector{Float64}
+    discharge_cost::Vector{Float64}
+    charge_efficiency::Vector{Float64}
+    discharge_efficiency::Vector{Float64}
+    loss_factor::Vector{Float64}
+    min_charge_rate::Vector{Float64}
+    max_charge_rate::Vector{Float64}
+    min_discharge_rate::Vector{Float64}
+    max_discharge_rate::Vector{Float64}
+    initial_level::Float64
+    min_ending_level::Float64
+    max_ending_level::Float64
+end
+
 Base.@kwdef mutable struct UnitCommitmentScenario
     buses_by_name::Dict{AbstractString,Bus}
     buses::Vector{Bus}
@@ -103,6 +124,8 @@ Base.@kwdef mutable struct UnitCommitmentScenario
     reserves::Vector{Reserve}
     thermal_units_by_name::Dict{AbstractString,ThermalUnit}
     thermal_units::Vector{ThermalUnit}
+    storage_units_by_name::Dict{AbstractString,StorageUnit}
+    storage_units::Vector{StorageUnit}
     time::Int
     time_step::Int
 end
diff --git a/src/model/build.jl b/src/model/build.jl
index a812011..0e27d42 100644
--- a/src/model/build.jl
+++ b/src/model/build.jl
@@ -99,6 +99,9 @@ function build_model(;
             for pu in sc.profiled_units
                 _add_profiled_unit!(model, pu, sc)
             end
+            for su in sc.storage_units
+                _add_storage_unit!(model, su, sc)
+            end
             _add_system_wide_eqs!(model, sc)
         end
         @objective(model, Min, model[:obj])
diff --git a/src/model/formulations/base/storage.jl b/src/model/formulations/base/storage.jl
new file mode 100644
index 0000000..93cf297
--- /dev/null
+++ b/src/model/formulations/base/storage.jl
@@ -0,0 +1,125 @@
+# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+function _add_storage_unit!(
+    model::JuMP.Model,
+    su::StorageUnit,
+    sc::UnitCommitmentScenario,
+)::Nothing
+    # Initialize variables
+    storage_level = _init(model, :storage_level)
+    charge_rate = _init(model, :charge_rate)
+    discharge_rate = _init(model, :discharge_rate)
+    is_charging = _init(model, :is_charging)
+    is_discharging = _init(model, :is_discharging)
+    eq_min_charge_rate = _init(model, :eq_min_charge_rate)
+    eq_max_charge_rate = _init(model, :eq_max_charge_rate)
+    eq_min_discharge_rate = _init(model, :eq_min_discharge_rate)
+    eq_max_discharge_rate = _init(model, :eq_max_discharge_rate)
+    # Initialize constraints
+    net_injection = _init(model, :expr_net_injection)
+    eq_storage_transition = _init(model, :eq_storage_transition)
+    eq_ending_level = _init(model, :eq_ending_level)
+    # time in hours
+    time_step = sc.time_step / 60
+
+    for t in 1:model[:instance].time
+        # Decision variable
+        storage_level[sc.name, su.name, t] = @variable(
+            model,
+            lower_bound = su.min_level[t],
+            upper_bound = su.max_level[t]
+        )
+        charge_rate[sc.name, su.name, t] = @variable(model)
+        discharge_rate[sc.name, su.name, t] = @variable(model)
+        is_charging[sc.name, su.name, t] = @variable(model, binary = true)
+        is_discharging[sc.name, su.name, t] = @variable(model, binary = true)
+
+        # Objective function terms ##### CHECK & FIXME
+        add_to_expression!(
+            model[:obj],
+            charge_rate[sc.name, su.name, t],
+            su.charge_cost[t] * sc.probability,
+        )
+
+        add_to_expression!(
+            model[:obj],
+            discharge_rate[sc.name, su.name, t],
+            su.discharge_cost[t] * sc.probability,
+        )
+
+        # Net injection
+        add_to_expression!(
+            net_injection[sc.name, su.bus.name, t],
+            discharge_rate[sc.name, su.name, t],
+            1.0,
+        )
+        add_to_expression!(
+            net_injection[sc.name, su.bus.name, t],
+            charge_rate[sc.name, su.name, t],
+            -1.0,
+        )
+
+        # Simultaneous charging and discharging
+        if !su.simultaneous_charge_and_discharge[t]
+            # Initialize the model dictionary
+            eq_simultaneous_charge_and_discharge =
+                _init(model, :eq_simultaneous_charge_and_discharge)
+            # Constraints
+            eq_simultaneous_charge_and_discharge[sc.name, su.name, t] =
+                @constraint(
+                    model,
+                    is_charging[sc.name, su.name, t] +
+                    is_discharging[sc.name, su.name, t] <= 1.0
+                )
+        end
+
+        # Charge and discharge constraints
+        eq_min_charge_rate[sc.name, su.name, t] = @constraint(
+            model,
+            charge_rate[sc.name, su.name, t] >=
+            is_charging[sc.name, su.name, t] * su.min_charge_rate[t]
+        )
+        eq_max_charge_rate[sc.name, su.name, t] = @constraint(
+            model,
+            charge_rate[sc.name, su.name, t] <=
+            is_charging[sc.name, su.name, t] * su.max_charge_rate[t]
+        )
+        eq_min_discharge_rate[sc.name, su.name, t] = @constraint(
+            model,
+            discharge_rate[sc.name, su.name, t] >=
+            is_discharging[sc.name, su.name, t] * su.min_discharge_rate[t]
+        )
+        eq_max_discharge_rate[sc.name, su.name, t] = @constraint(
+            model,
+            discharge_rate[sc.name, su.name, t] <=
+            is_discharging[sc.name, su.name, t] * su.max_discharge_rate[t]
+        )
+
+        # Storage energy transition constraint 
+        prev_storage_level =
+            t == 1 ? su.initial_level : storage_level[sc.name, su.name, t-1]
+        eq_storage_transition[sc.name, su.name, t] = @constraint(
+            model,
+            storage_level[sc.name, su.name, t] ==
+            (1 - su.loss_factor[t]) * prev_storage_level +
+            charge_rate[sc.name, su.name, t] *
+            time_step *
+            su.charge_efficiency[t] -
+            discharge_rate[sc.name, su.name, t] * time_step /
+            su.discharge_efficiency[t]
+        )
+
+        # Storage ending level constraint 
+        if t == sc.time
+            eq_ending_level[sc.name, su.name] = @constraint(
+                model,
+                su.min_ending_level <=
+                storage_level[sc.name, su.name, t] <=
+                su.max_ending_level
+            )
+        end
+    end
+    return
+end
diff --git a/src/solution/solution.jl b/src/solution/solution.jl
index 3cb6f41..c3a9cdf 100644
--- a/src/solution/solution.jl
+++ b/src/solution/solution.jl
@@ -103,6 +103,30 @@ function solution(model::JuMP.Model)::OrderedDict
                 ] for pu in sc.profiled_units
             )
         end
+        if !isempty(sc.storage_units)
+            sol[sc.name]["Storage level (MWh)"] =
+                timeseries(model[:storage_level], sc.storage_units, sc = sc)
+            sol[sc.name]["Is charging"] =
+                timeseries(model[:is_charging], sc.storage_units, sc = sc)
+            sol[sc.name]["Storage charging rates (MW)"] =
+                timeseries(model[:charge_rate], sc.storage_units, sc = sc)
+            sol[sc.name]["Storage charging cost (\$)"] = OrderedDict(
+                su.name => [
+                    value(model[:charge_rate][sc.name, su.name, t]) *
+                    su.charge_cost[t] for t in 1:instance.time
+                ] for su in sc.storage_units
+            )
+            sol[sc.name]["Is discharging"] =
+                timeseries(model[:is_discharging], sc.storage_units, sc = sc)
+            sol[sc.name]["Storage discharging rates (MW)"] =
+                timeseries(model[:discharge_rate], sc.storage_units, sc = sc)
+            sol[sc.name]["Storage discharging cost (\$)"] = OrderedDict(
+                su.name => [
+                    value(model[:discharge_rate][sc.name, su.name, t]) *
+                    su.discharge_cost[t] for t in 1:instance.time
+                ] for su in sc.storage_units
+            )
+        end
         sol[sc.name]["Spinning reserve (MW)"] = OrderedDict(
             r.name => OrderedDict(
                 g.name => [
diff --git a/src/transform/randomize/XavQiuAhm2021.jl b/src/transform/randomize/XavQiuAhm2021.jl
index f58d9cc..8404388 100644
--- a/src/transform/randomize/XavQiuAhm2021.jl
+++ b/src/transform/randomize/XavQiuAhm2021.jl
@@ -137,6 +137,11 @@ function _randomize_costs(
         α = rand(rng, distribution)
         pu.cost *= α
     end
+    for su in sc.storage_units
+        α = rand(rng, distribution)
+        su.charge_cost *= α
+        su.discharge_cost *= α
+    end
     return
 end
 
diff --git a/src/transform/slice.jl b/src/transform/slice.jl
index 3fbbdaa..ef90082 100644
--- a/src/transform/slice.jl
+++ b/src/transform/slice.jl
@@ -56,6 +56,21 @@ function slice(
             ps.demand = ps.demand[range]
             ps.revenue = ps.revenue[range]
         end
+        for su in sc.storage_units
+            su.min_level = su.min_level[range]
+            su.max_level = su.max_level[range]
+            su.simultaneous_charge_and_discharge =
+                su.simultaneous_charge_and_discharge[range]
+            su.charge_cost = su.charge_cost[range]
+            su.discharge_cost = su.discharge_cost[range]
+            su.charge_efficiency = su.charge_efficiency[range]
+            su.discharge_efficiency = su.discharge_efficiency[range]
+            su.loss_factor = su.loss_factor[range]
+            su.min_charge_rate = su.min_charge_rate[range]
+            su.max_charge_rate = su.max_charge_rate[range]
+            su.min_discharge_rate = su.min_discharge_rate[range]
+            su.max_discharge_rate = su.max_discharge_rate[range]
+        end
     end
     return modified
 end
diff --git a/src/validation/validate.jl b/src/validation/validate.jl
index 863b675..3310901 100644
--- a/src/validation/validate.jl
+++ b/src/validation/validate.jl
@@ -334,6 +334,195 @@ function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
                 end
             end
         end
+        for su in sc.storage_units
+            storage_level = solution[sc.name]["Storage level (MWh)"][su.name]
+            charge_rate =
+                solution[sc.name]["Storage charging rates (MW)"][su.name]
+            discharge_rate =
+                solution[sc.name]["Storage discharging rates (MW)"][su.name]
+            actual_charge_cost =
+                solution[sc.name]["Storage charging cost (\$)"][su.name]
+            actual_discharge_cost =
+                solution[sc.name]["Storage discharging cost (\$)"][su.name]
+            is_charging = bin(solution[sc.name]["Is charging"][su.name])
+            is_discharging = bin(solution[sc.name]["Is discharging"][su.name])
+            # time in hours
+            time_step = sc.time_step / 60
+
+            for t in 1:instance.time
+                # Unit must store at least its minimum level 
+                if storage_level[t] < su.min_level[t] - tol
+                    @error @sprintf(
+                        "Storage unit %s stores below its minimum level at time %d (%.2f < %.2f)",
+                        su.name,
+                        t,
+                        storage_level[t],
+                        su.min_level[t]
+                    )
+                    err_count += 1
+                end
+                # Unit must store at most its maximum level 
+                if storage_level[t] > su.max_level[t] + tol
+                    @error @sprintf(
+                        "Storage unit %s stores above its maximum level at time %d (%.2f > %.2f)",
+                        su.name,
+                        t,
+                        storage_level[t],
+                        su.max_level[t]
+                    )
+                    err_count += 1
+                end
+
+                if t == instance.time
+                    # Unit must store at least its minimum level at last time period
+                    if storage_level[t] < su.min_ending_level - tol
+                        @error @sprintf(
+                            "Storage unit %s stores below its minimum ending level (%.2f < %.2f)",
+                            su.name,
+                            storage_level[t],
+                            su.min_ending_level
+                        )
+                        err_count += 1
+                    end
+                    # Unit must store at most its maximum level at last time period
+                    if storage_level[t] > su.max_ending_level + tol
+                        @error @sprintf(
+                            "Storage unit %s stores above its maximum ending level (%.2f > %.2f)",
+                            su.name,
+                            storage_level[t],
+                            su.max_ending_level
+                        )
+                        err_count += 1
+                    end
+                end
+
+                # Unit must follow the energy transition constraint 
+                prev_level = t == 1 ? su.initial_level : storage_level[t-1]
+                current_level =
+                    (1 - su.loss_factor[t]) * prev_level +
+                    time_step * (
+                        charge_rate[t] * su.charge_efficiency[t] -
+                        discharge_rate[t] / su.discharge_efficiency[t]
+                    )
+                if abs(storage_level[t] - current_level) > tol
+                    @error @sprintf(
+                        "Storage unit %s has unexpected level at time %d (%.2f should be %.2f)",
+                        unit.name,
+                        t,
+                        storage_level[t],
+                        current_level
+                    )
+                    err_count += 1
+                end
+
+                # Unit cannot simultaneous charge and discharge if it is not allowed
+                if !su.simultaneous_charge_and_discharge[t] &&
+                   is_charging[t] &&
+                   is_discharging[t]
+                    @error @sprintf(
+                        "Storage unit %s is charging and discharging simultaneous at time %d",
+                        su.name,
+                        t
+                    )
+                    err_count += 1
+                end
+
+                # Unit must charge at least its minimum rate 
+                if is_charging[t] &&
+                   (charge_rate[t] < su.min_charge_rate[t] - tol)
+                    @error @sprintf(
+                        "Storage unit %s charges below its minimum limit at time %d (%.2f < %.2f)",
+                        unit.name,
+                        t,
+                        charge_rate[t],
+                        su.min_charge_rate[t]
+                    )
+                    err_count += 1
+                end
+                # Unit must charge at most its maximum rate 
+                if is_charging[t] &&
+                   (charge_rate[t] > su.max_charge_rate[t] + tol)
+                    @error @sprintf(
+                        "Storage unit %s charges above its maximum limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        charge_rate[t],
+                        su.max_charge_rate[t]
+                    )
+                    err_count += 1
+                end
+                # Unit must have zero charge when it is not charging
+                if !is_charging[t] && (charge_rate[t] > tol)
+                    @error @sprintf(
+                        "Storage unit %s charges power at time %d while not charging (%.2f > 0)",
+                        unit.name,
+                        t,
+                        charge_rate[t]
+                    )
+                    err_count += 1
+                end
+
+                # Unit must discharge at least its minimum rate 
+                if is_discharging[t] &&
+                   (discharge_rate[t] < su.min_discharge_rate[t] - tol)
+                    @error @sprintf(
+                        "Storage unit %s discharges below its minimum limit at time %d (%.2f < %.2f)",
+                        unit.name,
+                        t,
+                        discharge_rate[t],
+                        su.min_discharge_rate[t]
+                    )
+                    err_count += 1
+                end
+                # Unit must discharge at most its maximum rate 
+                if is_discharging[t] &&
+                   (discharge_rate[t] > su.max_discharge_rate[t] + tol)
+                    @error @sprintf(
+                        "Storage unit %s discharges above its maximum limit at time %d (%.2f > %.2f)",
+                        unit.name,
+                        t,
+                        discharge_rate[t],
+                        su.max_discharge_rate[t]
+                    )
+                    err_count += 1
+                end
+                # Unit must have zero discharge when it is not charging
+                if !is_discharging[t] && (discharge_rate[t] > tol)
+                    @error @sprintf(
+                        "Storage unit %s discharges power at time %d while not discharging (%.2f > 0)",
+                        unit.name,
+                        t,
+                        discharge_rate[t]
+                    )
+                    err_count += 1
+                end
+
+                # Compute storage costs 
+                charge_cost = su.charge_cost[t] * charge_rate[t]
+                discharge_cost = su.discharge_cost[t] * discharge_rate[t]
+                # Compare costs
+                if abs(actual_charge_cost[t] - charge_cost) > tol
+                    @error @sprintf(
+                        "Storage unit %s has unexpected charge cost at time %d (%.2f should be %.2f)",
+                        unit.name,
+                        t,
+                        actual_charge_cost[t],
+                        charge_cost
+                    )
+                    err_count += 1
+                end
+                if abs(actual_discharge_cost[t] - discharge_cost) > tol
+                    @error @sprintf(
+                        "Storage unit %s has unexpected discharge cost at time %d (%.2f should be %.2f)",
+                        unit.name,
+                        t,
+                        actual_discharge_cost[t],
+                        discharge_cost
+                    )
+                    err_count += 1
+                end
+            end
+        end
     end
     return err_count
 end
@@ -346,6 +535,8 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
             fixed_load = sum(b.load[t] for b in sc.buses)
             ps_load = 0
             production = 0
+            storage_charge = 0
+            storage_discharge = 0
             if length(sc.price_sensitive_loads) > 0
                 ps_load = sum(
                     solution[sc.name]["Price-sensitive loads (MW)"][ps.name][t]
@@ -364,23 +555,38 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
                     for pu in sc.profiled_units
                 )
             end
+            if length(sc.storage_units) > 0
+                storage_charge += sum(
+                    solution[sc.name]["Storage charging rates (MW)"][su.name][t]
+                    for su in sc.storage_units
+                )
+                storage_discharge += sum(
+                    solution[sc.name]["Storage discharging rates (MW)"][su.name][t]
+                    for su in sc.storage_units
+                )
+            end
             if "Load curtail (MW)" in keys(solution)
                 load_curtail = sum(
                     solution[sc.name]["Load curtail (MW)"][b.name][t] for
                     b in sc.buses
                 )
             end
-            balance = fixed_load - load_curtail - production + ps_load
+            balance =
+                fixed_load - load_curtail - production +
+                ps_load +
+                storage_charge - storage_discharge
 
             # Verify that production equals demand
             if abs(balance) > tol
                 @error @sprintf(
-                    "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)",
+                    "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f + %.2f - %.2f != 0)",
                     t,
                     fixed_load,
                     ps_load,
                     load_curtail,
                     production,
+                    storage_charge,
+                    storage_discharge,
                 )
                 err_count += 1
             end
diff --git a/test/fixtures/case14-storage.json.gz b/test/fixtures/case14-storage.json.gz
new file mode 100644
index 0000000..6a5451a
Binary files /dev/null and b/test/fixtures/case14-storage.json.gz differ
diff --git a/test/src/instance/read_test.jl b/test/src/instance/read_test.jl
index c9ae1fd..2ddffbc 100644
--- a/test/src/instance/read_test.jl
+++ b/test/src/instance/read_test.jl
@@ -139,20 +139,20 @@ function instance_read_test()
         sc = instance.scenarios[1]
         @test length(sc.profiled_units) == 2
 
-        first_pu = sc.profiled_units[1]
-        @test first_pu.name == "g7"
-        @test first_pu.bus.name == "b4"
-        @test first_pu.cost == [100.0 for t in 1:4]
-        @test first_pu.min_power == [60.0 for t in 1:4]
-        @test first_pu.max_power == [100.0 for t in 1:4]
+        pu1 = sc.profiled_units[1]
+        @test pu1.name == "g7"
+        @test pu1.bus.name == "b4"
+        @test pu1.cost == [100.0 for t in 1:4]
+        @test pu1.min_power == [60.0 for t in 1:4]
+        @test pu1.max_power == [100.0 for t in 1:4]
         @test sc.profiled_units_by_name["g7"].name == "g7"
 
-        second_pu = sc.profiled_units[2]
-        @test second_pu.name == "g8"
-        @test second_pu.bus.name == "b5"
-        @test second_pu.cost == [50.0 for t in 1:4]
-        @test second_pu.min_power == [0.0 for t in 1:4]
-        @test second_pu.max_power == [120.0 for t in 1:4]
+        pu2 = sc.profiled_units[2]
+        @test pu2.name == "g8"
+        @test pu2.bus.name == "b5"
+        @test pu2.cost == [50.0 for t in 1:4]
+        @test pu2.min_power == [0.0 for t in 1:4]
+        @test pu2.max_power == [120.0 for t in 1:4]
         @test sc.profiled_units_by_name["g8"].name == "g8"
     end
 
@@ -166,4 +166,64 @@ function instance_read_test()
         @test sc.thermal_units[6].commitment_status ==
               [false, nothing, true, nothing]
     end
+
+    @testset "read_benchmark storage" begin
+        instance = UnitCommitment.read(fixture("case14-storage.json.gz"))
+        sc = instance.scenarios[1]
+        @test length(sc.storage_units) == 4
+
+        su1 = sc.storage_units[1]
+        @test su1.name == "su1"
+        @test su1.bus.name == "b2"
+        @test su1.min_level == [0.0 for t in 1:4]
+        @test su1.max_level == [100.0 for t in 1:4]
+        @test su1.simultaneous_charge_and_discharge == [true for t in 1:4]
+        @test su1.charge_cost == [2.0 for t in 1:4]
+        @test su1.discharge_cost == [2.5 for t in 1:4]
+        @test su1.charge_efficiency == [1.0 for t in 1:4]
+        @test su1.discharge_efficiency == [1.0 for t in 1:4]
+        @test su1.loss_factor == [0.0 for t in 1:4]
+        @test su1.min_charge_rate == [0.0 for t in 1:4]
+        @test su1.max_charge_rate == [10.0 for t in 1:4]
+        @test su1.min_discharge_rate == [0.0 for t in 1:4]
+        @test su1.max_discharge_rate == [8.0 for t in 1:4]
+        @test su1.initial_level == 0.0
+        @test su1.min_ending_level == 0.0
+        @test su1.max_ending_level == 100.0
+        @test sc.storage_units_by_name["su1"].name == "su1"
+
+        su2 = sc.storage_units[2]
+        @test su2.name == "su2"
+        @test su2.bus.name == "b2"
+        @test su2.min_level == [10.0 for t in 1:4]
+        @test su2.simultaneous_charge_and_discharge == [false for t in 1:4]
+        @test su2.charge_cost == [3.0 for t in 1:4]
+        @test su2.discharge_cost == [3.5 for t in 1:4]
+        @test su2.charge_efficiency == [0.8 for t in 1:4]
+        @test su2.discharge_efficiency == [0.85 for t in 1:4]
+        @test su2.loss_factor == [0.01 for t in 1:4]
+        @test su2.min_charge_rate == [5.0 for t in 1:4]
+        @test su2.min_discharge_rate == [2.0 for t in 1:4]
+        @test su2.initial_level == 70.0
+        @test su2.min_ending_level == 80.0
+        @test su2.max_ending_level == 85.0
+        @test sc.storage_units_by_name["su2"].name == "su2"
+
+        su3 = sc.storage_units[3]
+        @test su3.bus.name == "b9"
+        @test su3.min_level == [10.0, 11.0, 12.0, 13.0]
+        @test su3.max_level == [100.0, 110.0, 120.0, 130.0]
+        @test su3.charge_cost == [2.0, 2.1, 2.2, 2.3]
+        @test su3.discharge_cost == [1.0, 1.1, 1.2, 1.3]
+        @test su3.charge_efficiency == [0.8, 0.81, 0.82, 0.82]
+        @test su3.discharge_efficiency == [0.85, 0.86, 0.87, 0.88]
+        @test su3.min_charge_rate == [5.0, 5.1, 5.2, 5.3]
+        @test su3.max_charge_rate == [10.0, 10.1, 10.2, 10.3]
+        @test su3.min_discharge_rate == [4.0, 4.1, 4.2, 4.3]
+        @test su3.max_discharge_rate == [8.0, 8.1, 8.2, 8.3]
+
+        su4 = sc.storage_units[4]
+        @test su4.simultaneous_charge_and_discharge ==
+              [false, false, true, true]
+    end
 end
diff --git a/test/src/transform/randomize/XavQiuAhm2021_test.jl b/test/src/transform/randomize/XavQiuAhm2021_test.jl
index 7819529..dfe397f 100644
--- a/test/src/transform/randomize/XavQiuAhm2021_test.jl
+++ b/test/src/transform/randomize/XavQiuAhm2021_test.jl
@@ -102,5 +102,28 @@ function transform_randomize_XavQiuAhm2021_test()
             test_approx(pu1.cost[1], 98.039)
             test_approx(pu2.cost[1], 48.385)
         end
+
+        @testset "storage unit cost" begin
+            sc = UnitCommitment.read(
+                fixture("case14-storage.json.gz"),
+            ).scenarios[1]
+            # Check original costs
+            su1 = sc.storage_units[1]
+            su3 = sc.storage_units[3]
+            test_approx(su1.charge_cost[4], 2.0)
+            test_approx(su1.discharge_cost[1], 2.5)
+            test_approx(su3.charge_cost[2], 2.1)
+            test_approx(su3.discharge_cost[3], 1.2)
+            randomize!(
+                sc,
+                XavQiuAhm2021.Randomization(randomize_load_profile = false),
+                rng = MersenneTwister(42),
+            )
+            # Check randomized costs
+            test_approx(su1.charge_cost[4], 1.961)
+            test_approx(su1.discharge_cost[1], 2.451)
+            test_approx(su3.charge_cost[2], 2.196)
+            test_approx(su3.discharge_cost[3], 1.255)
+        end
     end
 end
diff --git a/test/src/transform/slice_test.jl b/test/src/transform/slice_test.jl
index affe864..fdd86b4 100644
--- a/test/src/transform/slice_test.jl
+++ b/test/src/transform/slice_test.jl
@@ -65,4 +65,34 @@ function transform_slice_test()
             variable_names = true,
         )
     end
+
+    @testset "slice storage units" begin
+        instance = UnitCommitment.read(fixture("case14-storage.json.gz"))
+        modified = UnitCommitment.slice(instance, 2:4)
+        sc = modified.scenarios[1]
+
+        # Should update all time-dependent fields
+        for su in sc.storage_units
+            @test length(su.min_level) == 3
+            @test length(su.max_level) == 3
+            @test length(su.simultaneous_charge_and_discharge) == 3
+            @test length(su.charge_cost) == 3
+            @test length(su.discharge_cost) == 3
+            @test length(su.charge_efficiency) == 3
+            @test length(su.discharge_efficiency) == 3
+            @test length(su.loss_factor) == 3
+            @test length(su.min_charge_rate) == 3
+            @test length(su.max_charge_rate) == 3
+            @test length(su.min_discharge_rate) == 3
+            @test length(su.max_discharge_rate) == 3
+        end
+
+        # Should be able to build model without errors
+        optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
+        model = UnitCommitment.build_model(
+            instance = modified,
+            optimizer = optimizer,
+            variable_names = true,
+        )
+    end
 end