diff --git a/pywr-core/src/lib.rs b/pywr-core/src/lib.rs
index 02c79633..40e45056 100644
--- a/pywr-core/src/lib.rs
+++ b/pywr-core/src/lib.rs
@@ -79,6 +79,8 @@ pub enum PywrError {
     FlowConstraintsUndefined,
     #[error("storage constraints are undefined for this node")]
     StorageConstraintsUndefined,
+    #[error("can not add virtual storage node to a storage node")]
+    NoVirtualStorageOnStorageNode,
     #[error("timestep index out of range")]
     TimestepIndexOutOfRange,
     #[error("solver not initialised")]
diff --git a/pywr-core/src/model.rs b/pywr-core/src/model.rs
index 4e5ae3f6..fc7dd932 100644
--- a/pywr-core/src/model.rs
+++ b/pywr-core/src/model.rs
@@ -1352,12 +1352,13 @@ impl Model {
         min_volume: ConstraintValue,
         max_volume: ConstraintValue,
         reset: VirtualStorageReset,
+        cost: ConstraintValue,
     ) -> Result<VirtualStorageIndex, PywrError> {
         if let Ok(_agg_node) = self.get_virtual_storage_node_by_name(name, sub_name) {
             return Err(PywrError::NodeNameAlreadyExists(name.to_string()));
         }
 
-        let node_index = self.virtual_storage_nodes.push_new(
+        let vs_node_index = self.virtual_storage_nodes.push_new(
             name,
             sub_name,
             nodes,
@@ -1366,12 +1367,20 @@ impl Model {
             min_volume,
             max_volume,
             reset,
+            cost,
         );
 
+        // Link the virtual storage node to the nodes it is including
+        for node_idx in nodes {
+            let node = self.nodes.get_mut(node_idx)?;
+            node.add_virtual_storage(vs_node_index)?;
+        }
+
         // Add to the resolve order.
-        self.resolve_order.push(ComponentType::VirtualStorageNode(node_index));
+        self.resolve_order
+            .push(ComponentType::VirtualStorageNode(vs_node_index));
 
-        Ok(node_index)
+        Ok(vs_node_index)
     }
 
     /// Add a `parameters::Parameter` to the model
diff --git a/pywr-core/src/node.rs b/pywr-core/src/node.rs
index 17dbeac1..3cce421f 100644
--- a/pywr-core/src/node.rs
+++ b/pywr-core/src/node.rs
@@ -3,6 +3,7 @@ use crate::metric::Metric;
 use crate::model::Model;
 use crate::state::{NodeState, State};
 use crate::timestep::Timestep;
+use crate::virtual_storage::VirtualStorageIndex;
 use crate::PywrError;
 use std::ops::{Deref, DerefMut};
 
@@ -123,6 +124,13 @@ impl From<Metric> for ConstraintValue {
     }
 }
 
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum CostAggFunc {
+    Sum,
+    Max,
+    Min,
+}
+
 impl Node {
     /// Create a new input node
     pub fn new_input(node_index: &NodeIndex, name: &str, sub_name: Option<&str>) -> Self {
@@ -286,6 +294,24 @@ impl Node {
         }
     }
 
+    pub fn add_virtual_storage(&mut self, virtual_storage_index: VirtualStorageIndex) -> Result<(), PywrError> {
+        match self {
+            Self::Input(n) => {
+                n.cost.virtual_storage_nodes.push(virtual_storage_index);
+                Ok(())
+            }
+            Self::Output(n) => {
+                n.cost.virtual_storage_nodes.push(virtual_storage_index);
+                Ok(())
+            }
+            Self::Link(n) => {
+                n.cost.virtual_storage_nodes.push(virtual_storage_index);
+                Ok(())
+            }
+            Self::Storage(_) => Err(PywrError::NoVirtualStorageOnStorageNode),
+        }
+    }
+
     // /// Return a reference to a node's flow constraints if they exist.
     // fn flow_constraints(&self) -> Option<&FlowConstraints> {
     //     match self {
@@ -502,6 +528,17 @@ impl Node {
         }
     }
 
+    pub fn set_cost_agg_func(&mut self, agg_func: CostAggFunc) -> Result<(), PywrError> {
+        match self {
+            Self::Input(n) => n.set_cost_agg_func(agg_func),
+            Self::Link(n) => n.set_cost_agg_func(agg_func),
+            Self::Output(n) => n.set_cost_agg_func(agg_func),
+            Self::Storage(_) => return Err(PywrError::NoVirtualStorageOnStorageNode),
+        };
+
+        Ok(())
+    }
+
     pub fn get_outgoing_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
         match self {
             Self::Input(n) => n.get_cost(model, state),
@@ -628,10 +665,65 @@ impl StorageConstraints {
     }
 }
 
+/// Generic cost data for a node.
+#[derive(Debug, PartialEq)]
+struct NodeCost {
+    local: ConstraintValue,
+    virtual_storage_nodes: Vec<VirtualStorageIndex>,
+    agg_func: CostAggFunc,
+}
+
+impl Default for NodeCost {
+    fn default() -> Self {
+        Self {
+            local: ConstraintValue::None,
+            virtual_storage_nodes: Vec::new(),
+            agg_func: CostAggFunc::Max,
+        }
+    }
+}
+
+impl NodeCost {
+    fn get_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
+        let local_cost = match &self.local {
+            ConstraintValue::None => Ok(0.0),
+            ConstraintValue::Scalar(v) => Ok(*v),
+            ConstraintValue::Metric(m) => m.get_value(model, state),
+        }?;
+
+        let vs_costs: Vec<f64> = self
+            .virtual_storage_nodes
+            .iter()
+            .map(|idx| {
+                let vs = model.get_virtual_storage_node(idx)?;
+                vs.get_cost(model, state)
+            })
+            .collect::<Result<_, _>>()?;
+
+        let cost = match self.agg_func {
+            CostAggFunc::Sum => local_cost + vs_costs.iter().sum::<f64>(),
+            CostAggFunc::Max => local_cost.max(
+                vs_costs
+                    .into_iter()
+                    .max_by(|a, b| a.total_cmp(b))
+                    .unwrap_or(f64::NEG_INFINITY),
+            ),
+            CostAggFunc::Min => local_cost.min(
+                vs_costs
+                    .into_iter()
+                    .min_by(|a, b| a.total_cmp(b))
+                    .unwrap_or(f64::INFINITY),
+            ),
+        };
+
+        Ok(cost)
+    }
+}
+
 #[derive(Debug, PartialEq)]
 pub struct InputNode {
     pub meta: NodeMeta<NodeIndex>,
-    pub cost: ConstraintValue,
+    cost: NodeCost,
     pub flow_constraints: FlowConstraints,
     pub outgoing_edges: Vec<EdgeIndex>,
 }
@@ -640,20 +732,19 @@ impl InputNode {
     fn new(index: &NodeIndex, name: &str, sub_name: Option<&str>) -> Self {
         Self {
             meta: NodeMeta::new(index, name, sub_name),
-            cost: ConstraintValue::None,
+            cost: NodeCost::default(),
             flow_constraints: FlowConstraints::new(),
             outgoing_edges: Vec::new(),
         }
     }
     fn set_cost(&mut self, value: ConstraintValue) {
-        self.cost = value
+        self.cost.local = value
+    }
+    fn set_cost_agg_func(&mut self, agg_func: CostAggFunc) {
+        self.cost.agg_func = agg_func
     }
     fn get_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
-        match &self.cost {
-            ConstraintValue::None => Ok(0.0),
-            ConstraintValue::Scalar(v) => Ok(*v),
-            ConstraintValue::Metric(m) => m.get_value(model, state),
-        }
+        self.cost.get_cost(model, state)
     }
     fn set_min_flow(&mut self, value: ConstraintValue) {
         self.flow_constraints.min_flow = value;
@@ -678,7 +769,7 @@ impl InputNode {
 #[derive(Debug, PartialEq)]
 pub struct OutputNode {
     pub meta: NodeMeta<NodeIndex>,
-    pub cost: ConstraintValue,
+    cost: NodeCost,
     pub flow_constraints: FlowConstraints,
     pub incoming_edges: Vec<EdgeIndex>,
 }
@@ -687,20 +778,19 @@ impl OutputNode {
     fn new(index: &NodeIndex, name: &str, sub_name: Option<&str>) -> Self {
         Self {
             meta: NodeMeta::new(index, name, sub_name),
-            cost: ConstraintValue::None,
+            cost: NodeCost::default(),
             flow_constraints: FlowConstraints::new(),
             incoming_edges: Vec::new(),
         }
     }
     fn set_cost(&mut self, value: ConstraintValue) {
-        self.cost = value
+        self.cost.local = value
     }
     fn get_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
-        match &self.cost {
-            ConstraintValue::None => Ok(0.0),
-            ConstraintValue::Scalar(v) => Ok(*v),
-            ConstraintValue::Metric(m) => m.get_value(model, state),
-        }
+        self.cost.get_cost(model, state)
+    }
+    fn set_cost_agg_func(&mut self, agg_func: CostAggFunc) {
+        self.cost.agg_func = agg_func
     }
     fn set_min_flow(&mut self, value: ConstraintValue) {
         self.flow_constraints.min_flow = value;
@@ -725,7 +815,7 @@ impl OutputNode {
 #[derive(Debug, PartialEq)]
 pub struct LinkNode {
     pub meta: NodeMeta<NodeIndex>,
-    pub cost: ConstraintValue,
+    cost: NodeCost,
     pub flow_constraints: FlowConstraints,
     pub incoming_edges: Vec<EdgeIndex>,
     pub outgoing_edges: Vec<EdgeIndex>,
@@ -735,21 +825,20 @@ impl LinkNode {
     fn new(index: &NodeIndex, name: &str, sub_name: Option<&str>) -> Self {
         Self {
             meta: NodeMeta::new(index, name, sub_name),
-            cost: ConstraintValue::None,
+            cost: NodeCost::default(),
             flow_constraints: FlowConstraints::new(),
             incoming_edges: Vec::new(),
             outgoing_edges: Vec::new(),
         }
     }
     fn set_cost(&mut self, value: ConstraintValue) {
-        self.cost = value
+        self.cost.local = value
+    }
+    fn set_cost_agg_func(&mut self, agg_func: CostAggFunc) {
+        self.cost.agg_func = agg_func
     }
     fn get_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
-        match &self.cost {
-            ConstraintValue::None => Ok(0.0),
-            ConstraintValue::Scalar(v) => Ok(*v),
-            ConstraintValue::Metric(m) => m.get_value(model, state),
-        }
+        self.cost.get_cost(model, state)
     }
     fn set_min_flow(&mut self, value: ConstraintValue) {
         self.flow_constraints.min_flow = value;
diff --git a/pywr-core/src/virtual_storage.rs b/pywr-core/src/virtual_storage.rs
index 2c5c99f3..d16b5063 100644
--- a/pywr-core/src/virtual_storage.rs
+++ b/pywr-core/src/virtual_storage.rs
@@ -55,6 +55,7 @@ impl VirtualStorageVec {
         min_volume: ConstraintValue,
         max_volume: ConstraintValue,
         reset: VirtualStorageReset,
+        cost: ConstraintValue,
     ) -> VirtualStorageIndex {
         let node_index = VirtualStorageIndex(self.nodes.len());
         let node = VirtualStorage::new(
@@ -67,6 +68,7 @@ impl VirtualStorageVec {
             min_volume,
             max_volume,
             reset,
+            cost,
         );
         self.nodes.push(node);
         node_index
@@ -88,6 +90,7 @@ pub struct VirtualStorage {
     pub initial_volume: StorageInitialVolume,
     pub storage_constraints: StorageConstraints,
     pub reset: VirtualStorageReset,
+    pub cost: ConstraintValue,
 }
 
 impl VirtualStorage {
@@ -101,6 +104,7 @@ impl VirtualStorage {
         min_volume: ConstraintValue,
         max_volume: ConstraintValue,
         reset: VirtualStorageReset,
+        cost: ConstraintValue,
     ) -> Self {
         Self {
             meta: NodeMeta::new(index, name, sub_name),
@@ -110,6 +114,7 @@ impl VirtualStorage {
             initial_volume,
             storage_constraints: StorageConstraints::new(min_volume, max_volume),
             reset,
+            cost,
         }
     }
 
@@ -139,6 +144,14 @@ impl VirtualStorage {
         VirtualStorageState::new(0.0)
     }
 
+    pub fn get_cost(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
+        match &self.cost {
+            ConstraintValue::None => Ok(0.0),
+            ConstraintValue::Scalar(v) => Ok(*v),
+            ConstraintValue::Metric(m) => m.get_value(model, state),
+        }
+    }
+
     pub fn before(&self, timestep: &Timestep, model: &Model, state: &mut State) -> Result<(), PywrError> {
         let do_reset = if timestep.is_first() {
             // Set the initial volume if it is the first timestep.
@@ -219,12 +232,13 @@ mod tests {
     use crate::metric::Metric;
     use crate::model::Model;
     use crate::node::{ConstraintValue, StorageInitialVolume};
-    use crate::recorders::AssertionFnRecorder;
+    use crate::recorders::{AssertionFnRecorder, AssertionRecorder};
     use crate::scenario::ScenarioIndex;
     use crate::solvers::{ClpSolver, ClpSolverSettings};
-    use crate::test_utils::{default_timestepper, run_all_solvers};
+    use crate::test_utils::{default_timestepper, run_all_solvers, simple_model};
     use crate::timestep::Timestep;
     use crate::virtual_storage::{months_since_last_reset, VirtualStorageReset};
+    use ndarray::Array;
     use time::macros::date;
 
     /// Test the calculation of number of months since last reset
@@ -278,6 +292,7 @@ mod tests {
             ConstraintValue::Scalar(0.0),
             ConstraintValue::Scalar(100.0),
             VirtualStorageReset::Never,
+            ConstraintValue::Scalar(0.0),
         );
 
         // Setup a demand on output-0 and output-1
@@ -319,4 +334,35 @@ mod tests {
         // Test all solvers
         run_all_solvers(&model, &timestepper);
     }
+
+    #[test]
+    /// Test virtual storage node costs
+    fn test_virtual_storage_node_costs() {
+        let mut model = simple_model(1);
+        let timestepper = default_timestepper();
+
+        let nodes = vec![model.get_node_index_by_name("input", None).unwrap()];
+        // Virtual storage node cost is high enough to prevent any flow
+        model
+            .add_virtual_storage_node(
+                "vs",
+                None,
+                &nodes,
+                None,
+                StorageInitialVolume::Proportional(1.0),
+                ConstraintValue::Scalar(0.0),
+                ConstraintValue::Scalar(100.0),
+                VirtualStorageReset::Never,
+                ConstraintValue::Scalar(20.0),
+            )
+            .unwrap();
+
+        let expected = Array::zeros((366, 1));
+        let idx = model.get_node_by_name("output", None).unwrap().index();
+        let recorder = AssertionRecorder::new("output-flow", Metric::NodeInFlow(idx), expected, None, None);
+        model.add_recorder(Box::new(recorder)).unwrap();
+
+        // Test all solvers
+        run_all_solvers(&model, &timestepper);
+    }
 }
diff --git a/pywr-schema/src/nodes/annual_virtual_storage.rs b/pywr-schema/src/nodes/annual_virtual_storage.rs
index 5a054db9..3b11b08b 100644
--- a/pywr-schema/src/nodes/annual_virtual_storage.rs
+++ b/pywr-schema/src/nodes/annual_virtual_storage.rs
@@ -54,6 +54,11 @@ impl AnnualVirtualStorageNode {
             return Err(SchemaError::MissingInitialVolume(self.meta.name.to_string()));
         };
 
+        let cost = match &self.cost {
+            Some(v) => v.load(model, tables, data_path)?.into(),
+            None => ConstraintValue::Scalar(0.0),
+        };
+
         let min_volume = match &self.min_volume {
             Some(v) => v.load(model, tables, data_path)?.into(),
             None => ConstraintValue::Scalar(0.0),
@@ -84,6 +89,7 @@ impl AnnualVirtualStorageNode {
             min_volume,
             max_volume,
             reset,
+            cost,
         )?;
         Ok(())
     }
diff --git a/pywr-schema/src/nodes/monthly_virtual_storage.rs b/pywr-schema/src/nodes/monthly_virtual_storage.rs
index 0657372b..2c9bf47f 100644
--- a/pywr-schema/src/nodes/monthly_virtual_storage.rs
+++ b/pywr-schema/src/nodes/monthly_virtual_storage.rs
@@ -48,6 +48,11 @@ impl MonthlyVirtualStorageNode {
             return Err(SchemaError::MissingInitialVolume(self.meta.name.to_string()));
         };
 
+        let cost = match &self.cost {
+            Some(v) => v.load(model, tables, data_path)?.into(),
+            None => ConstraintValue::Scalar(0.0),
+        };
+
         let min_volume = match &self.min_volume {
             Some(v) => v.load(model, tables, data_path)?.into(),
             None => ConstraintValue::Scalar(0.0),
@@ -78,6 +83,7 @@ impl MonthlyVirtualStorageNode {
             min_volume,
             max_volume,
             reset,
+            cost,
         )?;
         Ok(())
     }
diff --git a/pywr-schema/src/nodes/virtual_storage.rs b/pywr-schema/src/nodes/virtual_storage.rs
index 671aadaf..30da54cb 100644
--- a/pywr-schema/src/nodes/virtual_storage.rs
+++ b/pywr-schema/src/nodes/virtual_storage.rs
@@ -36,6 +36,11 @@ impl VirtualStorageNode {
             return Err(SchemaError::MissingInitialVolume(self.meta.name.to_string()));
         };
 
+        let cost = match &self.cost {
+            Some(v) => v.load(model, tables, data_path)?.into(),
+            None => ConstraintValue::Scalar(0.0),
+        };
+
         let min_volume = match &self.min_volume {
             Some(v) => v.load(model, tables, data_path)?.into(),
             None => ConstraintValue::Scalar(0.0),
@@ -64,6 +69,7 @@ impl VirtualStorageNode {
             min_volume,
             max_volume,
             reset,
+            cost,
         )?;
         Ok(())
     }