diff --git a/temoa_model/temoa_initialize.py b/temoa_model/temoa_initialize.py
index ebb9c52e..750d8e05 100644
--- a/temoa_model/temoa_initialize.py
+++ b/temoa_model/temoa_initialize.py
@@ -1046,21 +1046,13 @@ def DemandActivityConstraintIndices ( M ):
 					yield (r,p,s,d,t,v,dem,s0,d0)
 
 def DemandConstraintIndices ( M ):
-	used_dems = set((r,dem) for r, p, dem in M.Demand.sparse_iterkeys())
-	DSD_keys = M.DemandSpecificDistribution.sparse_keys()
-	dem_slices = { (r,dem) : set(
-		(s, d)
-		for s in M.time_season
-		for d in M.time_of_day
-		if (r, s, d, dem) in DSD_keys )
-	  for (r,dem) in used_dems
-	}
 
 	indices = set(
 	  (r, p, s, d, dem)
 
 	  for r, p, dem in M.Demand.sparse_iterkeys()
-	  for s, d in dem_slices[ (r,dem) ]
+	  for s in M.time_season
+	  for d in M.time_of_day
 	)
 
 	return indices
diff --git a/temoa_model/temoa_model.py b/temoa_model/temoa_model.py
index 730c36be..15b6a30b 100755
--- a/temoa_model/temoa_model.py
+++ b/temoa_model/temoa_model.py
@@ -121,7 +121,7 @@ def temoa_create_model(name="Temoa"):
     # Define demand- and resource-related parameters
     M.DemandDefaultDistribution = Param(M.time_season, M.time_of_day, mutable=True)
     M.DemandSpecificDistribution = Param(
-        M.regions, M.time_season, M.time_of_day, M.commodity_demand, mutable=True
+        M.regions, M.time_season, M.time_of_day, M.commodity_demand, mutable=True, default=0
     )
 
     M.Demand = Param(M.regions, M.time_optimize, M.commodity_demand)
@@ -379,6 +379,10 @@ def temoa_create_model(name="Temoa"):
         M.StorageConstraints_rpsdtv, rule=StorageEnergy_Constraint
     )
 
+    # We make use of this following set in some of the storage constraints.
+    # Pre-computing it is considerably faster.
+    M.SegFracPerSeason = Param(M.time_season, initialize=SegFracPerSeason_rule)
+
     M.StorageEnergyUpperBoundConstraint = Constraint(
         M.StorageConstraints_rpsdtv, rule=StorageEnergyUpperBound_Constraint
     )
diff --git a/temoa_model/temoa_rules.py b/temoa_model/temoa_rules.py
index 5682f4c4..17a973bf 100644
--- a/temoa_model/temoa_rules.py
+++ b/temoa_model/temoa_rules.py
@@ -425,8 +425,6 @@ def Demand_Constraint(M, r, p, s, d, dem):
 could satisfy both an end-use and internal system demand, then the output from
 :math:`\textbf{FO}` and :math:`\textbf{FOA}` would be double counted.
 """
-    if (r,s,d,dem) not in M.DemandSpecificDistribution.sparse_keys():
-        return Constraint.Skip
 
     supply = sum(
         M.V_FlowOut[r, p, s, d, S_i, S_t, S_v, dem]
@@ -443,6 +441,7 @@ def Demand_Constraint(M, r, p, s, d, dem):
     DemandConstraintErrorCheck(supply + supply_annual, r, p, s, d, dem)
 
     expr = supply + supply_annual == M.Demand[r, p, dem] * M.DemandSpecificDistribution[r, s, d, dem]
+
     return expr
 
 def DemandActivity_Constraint(M, r, p, s, d, t, v, dem, s_0, d_0):
@@ -474,9 +473,6 @@ def DemandActivity_Constraint(M, r, p, s, d, t, v, dem, s_0, d_0):
 variations, and therefore the equation above only includes :math:`\textbf{FO}`
 and not  :math:`\textbf{FOA}`
 """
-    if (r,s,d,dem) not in M.DemandSpecificDistribution.sparse_keys():
-        return Constraint.Skip
-    DSD = M.DemandSpecificDistribution  # lazy programmer
 
     act_a = sum(
         M.V_FlowOut[r, p, s_0, d_0, S_i, t, v, dem]
@@ -487,7 +483,7 @@ def DemandActivity_Constraint(M, r, p, s, d, t, v, dem, s_0, d_0):
         for S_i in M.ProcessInputsByOutput[r, p, t, v, dem]
     )
 
-    expr = act_a * DSD[r, s, d, dem] == act_b * DSD[r, s_0, d_0, dem]
+    expr = act_a * M.DemandSpecificDistribution[r, s, d, dem] == act_b * M.DemandSpecificDistribution[r, s_0, d_0, dem]
     return expr
 
 
@@ -989,7 +985,7 @@ def StorageEnergyUpperBound_Constraint(M, r, p, s, d, t, v):
         M.V_Capacity[r, t, v]
         * M.CapacityToActivity[r, t]
         * (M.StorageDuration[r, t] / 8760)
-        * sum(M.SegFrac[s,S_d] for S_d in M.time_of_day) * 365
+        * M.SegFracPerSeason[s] * 365
         * value(M.ProcessLifeFrac[r, p, t, v])
     )
     expr = M.V_StorageLevel[r, p, s, d, t, v] <= energy_capacity
@@ -2600,7 +2596,8 @@ def ParamLoanAnnualize_rule(M, r, t, v):
 
     return annualized_rate
 
-
+def SegFracPerSeason_rule(M, s):
+    return sum(M.SegFrac[s, S_d] for S_d in M.time_of_day)
 
 def LinkedEmissionsTech_Constraint(M, r, p, s, d, t, v, e):
     r"""
diff --git a/temoa_model/temoa_run.py b/temoa_model/temoa_run.py
index f6291a92..a81aea82 100755
--- a/temoa_model/temoa_run.py
+++ b/temoa_model/temoa_run.py
@@ -378,7 +378,23 @@ def solve_temoa_instance (self):
 				if self.options.neos:
 				    self.result = self.optimizer.solve(self.instance, opt=self.options.solver)
 				else:
-				    self.result = self.optimizer.solve( self.instance, suffixes=['dual'],# 'rc', 'slack'],
+					if self.options.solver == 'cbc':
+						# Solver options. Reference: https://genxproject.github.io/GenX/dev/solver_configuration/
+						self.optimizer.options["dualTolerance"] = 1e-6
+						self.optimizer.options["primalTolerance"] = 1e-6
+						self.optimizer.options["zeroTolerance"] = 1e-12
+						self.optimizer.options["crossover"] = 'off'
+
+					elif self.options.solver == 'cplex':
+						# Note: these parameter values are taken to be the same as those in PyPSA (see: https://pypsa-eur.readthedocs.io/en/latest/configuration.html)
+						self.optimizer.options["lpmethod"] = 4 # barrier
+						self.optimizer.options["solutiontype"] = 2 # non basic solution, ie no crossover
+						self.optimizer.options["barrier convergetol"] = 1.e-5
+						self.optimizer.options["feasopt tolerance"] = 1.e-6
+
+					# Note: still need to add gurobi parameters.
+
+					self.result = self.optimizer.solve( self.instance, suffixes=['dual'],# 'rc', 'slack'],
 								keepfiles=self.options.keepPyomoLP,
 								symbolic_solver_labels=self.options.keepPyomoLP )
 				yield '\t\t\t\t\t\t[%8.2f]\n' % duration()