diff --git a/pareto_analysis/pareto_analysis.py b/pareto_analysis/pareto_analysis.py
index 2c93c4d9049a4722a8106b2ea0f34b5c3b8a2f2a..3a8f51b5d3784703491f334332f20df20785f472 100644
--- a/pareto_analysis/pareto_analysis.py
+++ b/pareto_analysis/pareto_analysis.py
@@ -26,7 +26,7 @@ from pyomo.opt import SolverStatus, TerminationCondition
 import matplotlib.pyplot as plot
 import os
 
-def pareto_analysis(prosumer, strategy_names):
+def pareto_analysis(prosumer, model, strategy_names):
     strategy_name_1 = strategy_names[0]
     strategy_name_2 = strategy_names[1]
     solver = pyomo.SolverFactory("gurobi")
@@ -34,50 +34,50 @@ def pareto_analysis(prosumer, strategy_names):
     solver.options['Presolve'] = 2
     solver.options['TimeLimit'] = 200
     payoff = [[0.0, 0.0], [0.0, 0.0]]
-    f1 = getattr(prosumer._model_sizing, 'f' + strategy_name_1)
-    O1 = getattr(prosumer._model_sizing, 'O' + strategy_name_1)
-    f2 = getattr(prosumer._model_sizing, 'f' + strategy_name_2)
-    O2 = getattr(prosumer._model_sizing, 'O' + strategy_name_2)
+    f1 = getattr(model, 'f' + strategy_name_1)
+    O1 = getattr(model, 'O' + strategy_name_1)
+    f2 = getattr(model, 'f' + strategy_name_2)
+    O2 = getattr(model, 'O' + strategy_name_2)
     O2.deactivate()
-    solver.solve(prosumer._model_sizing)
+    solver.solve(model)
     payoff[0][0] = pyomo.value(f1)
-    prosumer._model_sizing.cp11 = pyomo.Constraint(expr = f1 == payoff[0][0])
+    model.cp11 = pyomo.Constraint(expr = f1 == payoff[0][0])
     O1.deactivate()
     O2.activate()
-    solver.solve(prosumer._model_sizing)
+    solver.solve(model)
     payoff[0][1] = pyomo.value(f2)
-    prosumer._model_sizing.del_component(prosumer._model_sizing.cp11)
-    solver.solve(prosumer._model_sizing)
+    model.del_component(model.cp11)
+    solver.solve(model)
     payoff[1][1] = pyomo.value(f2)
-    prosumer._model_sizing.cp22 = pyomo.Constraint(expr = f2 == payoff[1][1])
+    model.cp22 = pyomo.Constraint(expr = f2 == payoff[1][1])
     O2.deactivate()
     O1.activate()
-    solver.solve(prosumer._model_sizing)
+    solver.solve(model)
     payoff[1][0] = pyomo.value(f1)
     r2 = payoff[1][1] - payoff[0][1]
     g2 = 4
     i2 = 0
     number_of_solutions = 0
-    prosumer._model_sizing.s2 = pyomo.Var(bounds = (0, None))
-    prosumer._model_sizing.e2 = pyomo.Param(mutable = True)
-    prosumer._model_sizing.del_component(prosumer._model_sizing.cp22)
+    model.s2 = pyomo.Var(bounds = (0, None))
+    model.e2 = pyomo.Param(mutable = True)
+    model.del_component(model.cp22)
     if O2.sense == pyomo.maximize:
         factor = -1
     else:
         factor = 1
-    prosumer._model_sizing.cp2 = pyomo.Constraint(expr = f2 + factor * prosumer._model_sizing.s2 == prosumer._model_sizing.e2)
+    model.cp2 = pyomo.Constraint(expr = f2 + factor * model.s2 == model.e2)
     O1.deactivate()
-    prosumer._model_sizing.O = pyomo.Objective(expr = f1 + 0.001 * (prosumer._model_sizing.s2 / r2), sense = O1.sense)
+    model.O = pyomo.Objective(expr = f1 + 0.001 * (model.s2 / r2), sense = O1.sense)
     f1_list = []
     f2_list = []
     pareto_rsl = dict()
     while i2 <= g2:
-        prosumer._model_sizing.e2 = payoff[0][1] + (i2 / g2) * r2
-        solver_result = solver.solve(prosumer._model_sizing)
+        model.e2 = payoff[0][1] + (i2 / g2) * r2
+        solver_result = solver.solve(model)
         if solver_result.solver.status == SolverStatus.ok and solver_result.solver.termination_condition == TerminationCondition.optimal:
             f1_list.append(pyomo.value(f1))
             f2_list.append(pyomo.value(f2))
-            pareto_rsl[number_of_solutions] = (pyomo.value(prosumer._model_sizing.x1), pyomo.value(prosumer._model_sizing.x2))
+            pareto_rsl[number_of_solutions] = (pyomo.value(model.x1), pyomo.value(model.x2))
             number_of_solutions += 1
             i2 += 1
         else: