diff --git a/.gitignore b/.gitignore
index f7a5d2f..5c84ab1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,33 @@
-## GitIgnore for HyPaD
-### MATLAB temporary files
-*.aux
-
-### Gurobi problem file 
-*.lp
\ No newline at end of file
+## GitIgnore for HyPaD
+.DS_Store
+
+### MATLAB temporary files
+*.aux
+
+### Gurobi problem file 
+*.lp
+
+### Python ignore
+# Byte-compiled / optimized / DLL files
+python/__pycache__/
+python/problems/__pycache__/
+python/solver/__pycache__/
+
+# PyInstaller
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# Spyder project settings
+.spyderproject
+.spyproject
\ No newline at end of file
diff --git a/README.md b/README.md
index c0c49fd..78fd441 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,10 @@
 # HyPaD
-This repository contains a MATLAB implementation of the Hybrid Patch Decomposition Algorithm to solve multi-objective mixed-integer convex optimization problems.
+This repository contains a Matlab and a Python implementation of the Hybrid Patch Decomposition Algorithm to solve multi-objective mixed-integer convex optimization problems.
 
 ## Getting Started
 You can start using HyPaD by downloading or cloning this repository using the green button near the top of the [GitHub page](https://github.com/LeoWarnow/HyPaD).
-Then, just open the [UserFile.m](https://github.com/LeoWarnow/HyPaD/blob/main/UserFile.m).
-It serves as an interface to the HyPaD algorithm and provides all instructions that you need to get started.
+Then, just open the [UserFile.m](https://github.com/LeoWarnow/HyPaD/blob/main/matlab/UserFile.m) to access the Matlab implementation or the [UserFile.ipynb](https://github.com/LeoWarnow/HyPaD/blob/main/python/UserFile.ipynb) for a Jupyter Notebook to access the Python implementation.
+Both serve as an interface to AdEnA and provide all instructions that you need to get started.
 
 ## References
 As this implementation is based on the scientific work by [Gabriele Eichfelder](https://www.tu-ilmenau.de/mmor/team/gabriele-eichfelder/) and [Leo Warnow](https://www.tu-ilmenau.de/mmor/team/leo-warnow/), please cite the corresponding paper when using this code:
diff --git a/UserFile.m b/matlab/UserFile.m
similarity index 100%
rename from UserFile.m
rename to matlab/UserFile.m
diff --git a/callSolver.m b/matlab/callSolver.m
similarity index 100%
rename from callSolver.m
rename to matlab/callSolver.m
diff --git a/initSession.m b/matlab/initSession.m
similarity index 100%
rename from initSession.m
rename to matlab/initSession.m
diff --git a/problems/H1.m b/matlab/problems/H1.m
similarity index 100%
rename from problems/H1.m
rename to matlab/problems/H1.m
diff --git a/problems/T1.m b/matlab/problems/T1.m
similarity index 100%
rename from problems/T1.m
rename to matlab/problems/T1.m
diff --git a/problems/T10.m b/matlab/problems/T10.m
similarity index 100%
rename from problems/T10.m
rename to matlab/problems/T10.m
diff --git a/problems/T2.m b/matlab/problems/T2.m
similarity index 100%
rename from problems/T2.m
rename to matlab/problems/T2.m
diff --git a/problems/T3.m b/matlab/problems/T3.m
similarity index 100%
rename from problems/T3.m
rename to matlab/problems/T3.m
diff --git a/problems/T4.m b/matlab/problems/T4.m
similarity index 100%
rename from problems/T4.m
rename to matlab/problems/T4.m
diff --git a/problems/T5.m b/matlab/problems/T5.m
similarity index 100%
rename from problems/T5.m
rename to matlab/problems/T5.m
diff --git a/problems/T6.m b/matlab/problems/T6.m
similarity index 100%
rename from problems/T6.m
rename to matlab/problems/T6.m
diff --git a/problems/T8.m b/matlab/problems/T8.m
similarity index 100%
rename from problems/T8.m
rename to matlab/problems/T8.m
diff --git a/problems/T9.m b/matlab/problems/T9.m
similarity index 100%
rename from problems/T9.m
rename to matlab/problems/T9.m
diff --git a/solver/FeasibilityCheckPatch.m b/matlab/solver/FeasibilityCheckPatch.m
similarity index 100%
rename from solver/FeasibilityCheckPatch.m
rename to matlab/solver/FeasibilityCheckPatch.m
diff --git a/solver/HyPaD.m b/matlab/solver/HyPaD.m
similarity index 100%
rename from solver/HyPaD.m
rename to matlab/solver/HyPaD.m
diff --git a/solver/PS.m b/matlab/solver/PS.m
similarity index 100%
rename from solver/PS.m
rename to matlab/solver/PS.m
diff --git a/solver/RSUP.m b/matlab/solver/RSUP.m
similarity index 100%
rename from solver/RSUP.m
rename to matlab/solver/RSUP.m
diff --git a/solver/SNIABoxesDynamic.m b/matlab/solver/SNIABoxesDynamic.m
similarity index 100%
rename from solver/SNIABoxesDynamic.m
rename to matlab/solver/SNIABoxesDynamic.m
diff --git a/solver/SNIABoxesFixed.m b/matlab/solver/SNIABoxesFixed.m
similarity index 100%
rename from solver/SNIABoxesFixed.m
rename to matlab/solver/SNIABoxesFixed.m
diff --git a/solver/SNIABoxesInit.m b/matlab/solver/SNIABoxesInit.m
similarity index 100%
rename from solver/SNIABoxesInit.m
rename to matlab/solver/SNIABoxesInit.m
diff --git a/solver/SNIAList.m b/matlab/solver/SNIAList.m
similarity index 100%
rename from solver/SNIAList.m
rename to matlab/solver/SNIAList.m
diff --git a/solver/SNIAListInit.m b/matlab/solver/SNIAListInit.m
similarity index 100%
rename from solver/SNIAListInit.m
rename to matlab/solver/SNIAListInit.m
diff --git a/solver/SUP.m b/matlab/solver/SUP.m
similarity index 100%
rename from solver/SUP.m
rename to matlab/solver/SUP.m
diff --git a/solver/WSNLP.m b/matlab/solver/WSNLP.m
similarity index 100%
rename from solver/WSNLP.m
rename to matlab/solver/WSNLP.m
diff --git a/solver/computeNDS.m b/matlab/solver/computeNDS.m
similarity index 100%
rename from solver/computeNDS.m
rename to matlab/solver/computeNDS.m
diff --git a/solver/plotBoxes.m b/matlab/solver/plotBoxes.m
similarity index 100%
rename from solver/plotBoxes.m
rename to matlab/solver/plotBoxes.m
diff --git a/solver/updateLLB3.m b/matlab/solver/updateLLB3.m
similarity index 100%
rename from solver/updateLLB3.m
rename to matlab/solver/updateLLB3.m
diff --git a/solver/updateLUB3.m b/matlab/solver/updateLUB3.m
similarity index 100%
rename from solver/updateLUB3.m
rename to matlab/solver/updateLUB3.m
diff --git a/problems/P1.m b/problems/P1.m
deleted file mode 100644
index 339f865..0000000
--- a/problems/P1.m
+++ /dev/null
@@ -1,52 +0,0 @@
-function [n,m,p,q,f,g,Df,Dg,Aineq,bineq,Aeq,beq,lb,ub,x0,is_convex,is_quadratic] = P1(~)
-%P1 A nonconvex quadratic test instance
-%   This example was taken from:
-%   Kristo Mela, Juhani Koski and Risto Silvennoinen. Algorithm for
-%   Generating the Pareto Optimal Set of Multiobjective Nonlinear
-%   Mixed-Integer Optimization Problems, AIAA (2007)
-
-% Dimension of decision and criterion space
-n = 2; % Continuous variables
-m = 8; % Integer variables
-p = 2; % Dimension criterion space
-q = 1; % Number of constraints
-
-% Problem type
-is_convex = false;
-is_quadratic = true;
-
-% Objective function
-G = [1,-1,2,0,0,0,0,0,0,0;...
-    -1,2,0,0,2,0,0,0,0,0;...
-    0,0,3,0,2,0,0,0,0,0;...
-    2,0,0,4,0,2,0,2,0,0;...
-    0,0,0,0,5,2,0,0,0,0;...
-    0,0,0,0,0,6,0,0,0,0;...
-    0,0,0,0,0,0,7,0,0,0;...
-    0,0,0,0,0,0,0,0,0,0;...
-    0,0,0,0,0,2,0,0,0,0;...
-    0,2,0,0,0,0,0,0,0,10];
-c1 = [-1;-1;1;-10;0;1;-2;0;3;0];
-c2 = [1;2;-1;1;5;-2;0;6;0;-3];
-f = @(x) [0.5*x'*G*x+c1'*x;c2'*x];
-Df = @(x) [(0.5*(G+G')*x+c1)';c2'];
-
-% Linear constraints (Aineq*x <= bineq, Aeq*x = beq)
-Aineq = [];
-bineq = [];
-Aeq = [];
-beq = [];
-
-% Lower and upper bounds (lb <= x <= ub)
-lb = [-1;-1;0;0;0;0;0;0;0;0];
-ub = ones(n+m,1);
-% z = [-10;-10];
-% Z = [25;15];
-
-% Start point x0
-x0 = ceil((lb+ub)/2);
-
-% Non-linear constraints (g(x) <= 0)
-g = @(x) 0;
-Dg = @(x) zeros(1,n+m);
-end
\ No newline at end of file
diff --git a/problems/P5.m b/problems/P5.m
deleted file mode 100644
index ed50337..0000000
--- a/problems/P5.m
+++ /dev/null
@@ -1,35 +0,0 @@
-function [n,m,p,q,f,g,Df,Dg,Aineq,bineq,Aeq,beq,lb,ub,x0,is_convex,is_quadratic] = P5(param)
-%P5 A simple scalable convex example
-
-% Dimension of decision and criterion space
-n = param(1); % Continuous variables
-m = n; % Integer variables
-p = n; % Dimension criterion space
-q = 2; % Number of constraints
-
-% Problem type
-is_convex = true;
-is_quadratic = true;
-
-% Objective function
-f = @(x) [x(1:n)+x(n+1:n+m)];
-Df = @(x) [eye(n),eye(m)];
-
-
-% Linear constraints (Aineq*x <= bineq, Aeq*x = beq)
-Aineq = [zeros(1,n),ones(1,m)];
-bineq = [3];
-Aeq = [];
-beq = [];
-
-% Lower and upper bounds (lb <= x <= ub)
-lb = -3.*ones(n+m,1);
-ub = 3.*ones(n+m,1);
-
-% Start point x0
-x0 = ceil((lb+ub)/2);
-
-% Non-linear constraints (g(x) <= 0)
-g = @(x) [sum(x(1:n).^2)-1;sum(x(n+1:n+m).^2)-11];
-Dg = @(x) [2.*x(1:n,1)',zeros(1,m);zeros(1,n),2.*x(n+1:n+m,1)'];
-end
\ No newline at end of file
diff --git a/python/UserFile.ipynb b/python/UserFile.ipynb
new file mode 100644
index 0000000..8467f2d
--- /dev/null
+++ b/python/UserFile.ipynb
@@ -0,0 +1,86 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# General packages\n",
+    "%matplotlib widget\n",
+    "import numpy as np\n",
+    "\n",
+    "# AdEnA\n",
+    "from solver.HyPaD import *\n",
+    "from solver.plot_functions import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Problem\n",
+    "from problems.H1 import *\n",
+    "params = [2,2]\n",
+    "\n",
+    "# Set quality epsilon and offset\n",
+    "EPSILON = 0.1\n",
+    "OFFSET = EPSILON*1e-3\n",
+    "\n",
+    "# Should the result be plotted (m = 2 and m = 3 only) [1 == yes, 0 == no]\n",
+    "plot_result = 1\n",
+    "\n",
+    "# Load optimization problem\n",
+    "model, model_parameters = build_model(params)\n",
+    "\n",
+    "# Initialize enclosure\n",
+    "L, U = init_enclosure(model, model_parameters, OFFSET)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Solve\n",
+    "L,U,N,ids,it,flag = HyPaD(build_model, params, L, U, EPSILON, OFFSET)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot\n",
+    "if plot_result == 1:\n",
+    "    plot_bounds(L,U,model_parameters.p)\n",
+    "    plot_boxes(L,U,model_parameters.p)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/python/problems/H1.py b/python/problems/H1.py
new file mode 100644
index 0000000..64a8026
--- /dev/null
+++ b/python/problems/H1.py
@@ -0,0 +1,82 @@
+"""
+H1 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(params):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = params[0] # Continuous variables
+    m = params[1] # Integer variables
+    p = 2         # Dimension criterion space
+    q = 1         # Number of constraints
+    assert n % 2 == 0, 'Number of continuous variables has to be even.'
+    assert m % 2 == 0, 'Number of integer variables has to be even.'
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-2,2)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = sum(model.x[i] for i in range(n // 2)) + sum(model.x[i] ** 2 for i in range(n, n + m // 2)) - sum(model.x[i] for i in range(n + m // 2, n + m))).deactivate()
+    model.f.add(expr = sum(model.x[i] for i in range(n // 2, n)) - sum(model.x[i] for i in range(n, n + m // 2)) + sum(model.x[i] ** 2 for i in range(n + m // 2, n + m))).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    for j in range(n//2):
+        model.Df[0, j] = 1
+        model.Df[1, j] = 0
+    for j in range(n//2, n):
+        model.Df[0, j] = 0
+        model.Df[1, j] = 1
+    for j in range(n,n+m//2):
+        model.Df[0, j] = 2 * model.x[j]
+        model.Df[1, j] = -1
+    for j in range(n+m//2,n+m):
+        model.Df[0, j] = -1
+        model.Df[1, j] = 2 * model.x[j]
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = sum(model.x[i] ** 2 for i in range(n)) - 1 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    for j in range(n):
+        model.Dg[0, j] = 2 * model.x[j]
+    for j in range(n, n+m):
+        model.Dg[0, j] = 0
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.empty((1,p))
+    model_parameters.U = np.empty((1,p))
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T1.py b/python/problems/T1.py
new file mode 100644
index 0000000..a02da7f
--- /dev/null
+++ b/python/problems/T1.py
@@ -0,0 +1,75 @@
+"""
+T1 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+
+It is taken from:
+Marianna De Santis, Gabriele Eichfelder, Julia Niebling, Stefan
+Rocktäschel. Solving Multiobjective Mixed Integer Convex Optimization
+Problems, SIAM Journal on Optimization (2020)
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 1         # Continuous variables
+    m = 1         # Integer variables
+    p = 2         # Dimension criterion space
+    q = 1         # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-4,4)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[1]).deactivate()
+    model.f.add(expr = model.x[0] ** 2 + model.x[1] ** 2).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 1
+    model.Df[1, 0] = 2 * model.x[0]  
+    model.Df[1, 1] = 2 * model.x[1]
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = (model.x[0] - 2) ** 2 + (model.x[1] - 2) ** 2 - 36 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    model.Dg[0, 0] = 2 * (model.x[0] - 2)
+    model.Dg[0, 1] = 2 * (model.x[1] - 2)
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-6,0]])
+    model_parameters.U = np.array([[6,20]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T10.py b/python/problems/T10.py
new file mode 100644
index 0000000..da05ceb
--- /dev/null
+++ b/python/problems/T10.py
@@ -0,0 +1,119 @@
+"""
+T10 is a bi-objective mixed-integer convex test instance with a
+non-quadratic objective function and quadratic constraint functions.
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 4         # Continuous variables
+    m = 4         # Integer variables
+    p = 2         # Dimension criterion space
+    q = 4         # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-20,20)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-20,20)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[2] + model.x[4] + exp(model.x[6]) - 1).deactivate()
+    model.f.add(expr = model.x[1] + model.x[3] + model.x[5] + model.x[7]).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 0
+    model.Df[0, 2] = 1
+    model.Df[0, 3] = 0
+    model.Df[0, 4] = 1
+    model.Df[0, 5] = 0
+    model.Df[0, 6] = exp(model.x[6])
+    model.Df[0, 7] = 0
+
+    model.Df[1, 0] = 0
+    model.Df[1, 1] = 1
+    model.Df[1, 2] = 0
+    model.Df[1, 3] = 1
+    model.Df[1, 4] = 0
+    model.Df[1, 5] = 1
+    model.Df[1, 6] = 0
+    model.Df[1, 7] = 1
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = model.x[0] ** 2 + model.x[1] ** 2 - 1 <= 0)
+    model.g.add(expr = model.x[2] ** 2 + model.x[3] ** 2 - 1 <= 0)
+    model.g.add(expr = (model.x[4] - 2) ** 2 + (model.x[5] - 5) ** 2 - 10 <= 0)
+    model.g.add(expr = (model.x[6] - 3) ** 2 + (model.x[7] - 8) ** 2 - 10 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    model.Dg[0, 0] = 2 * model.x[0]
+    model.Dg[0, 1] = 2 * model.x[1]
+    model.Dg[0, 2] = 0
+    model.Dg[0, 3] = 0
+    model.Dg[0, 4] = 0
+    model.Dg[0, 5] = 0
+    model.Dg[0, 6] = 0
+    model.Dg[0, 7] = 0
+
+    model.Dg[1, 0] = 0
+    model.Dg[1, 1] = 0
+    model.Dg[1, 2] = 2 * model.x[2]
+    model.Dg[1, 3] = 2 * model.x[3]
+    model.Dg[1, 4] = 0
+    model.Dg[1, 5] = 0
+    model.Dg[1, 6] = 0
+    model.Dg[1, 7] = 0
+
+    model.Dg[2, 0] = 0
+    model.Dg[2, 1] = 0
+    model.Dg[2, 2] = 0
+    model.Dg[2, 3] = 0
+    model.Dg[2, 4] = 2 * (model.x[4] - 2)
+    model.Dg[2, 5] = 2 * (model.x[5] - 5)
+    model.Dg[2, 6] = 0
+    model.Dg[2, 7] = 0
+
+    model.Dg[3, 0] = 0
+    model.Dg[3, 1] = 0
+    model.Dg[3, 2] = 0
+    model.Dg[3, 3] = 0
+    model.Dg[3, 4] = 0
+    model.Dg[3, 5] = 0
+    model.Dg[3, 6] = 2 * (model.x[6] - 3)
+    model.Dg[3, 7] = 2 * (model.x[7] - 8)
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-3,5]])
+    model_parameters.U = np.array([[12,22]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T3.py b/python/problems/T3.py
new file mode 100644
index 0000000..05ca5d5
--- /dev/null
+++ b/python/problems/T3.py
@@ -0,0 +1,79 @@
+"""
+T3 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+
+It is taken from:
+Marianna De Santis, Gabriele Eichfelder, Julia Niebling, Stefan
+Rocktäschel. Solving Multiobjective Mixed Integer Convex Optimization
+Problems, SIAM Journal on Optimization (2020)
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(params):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = params[0] # Continuous variables
+    m = params[1] # Integer variables
+    p = 2         # Dimension criterion space
+    q = 1         # Number of constraints
+    assert n % 2 == 0, 'Number of continuous variables has to be even.'
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-2,2)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0]).deactivate()
+    model.f.add(expr = model.x[1] + sum(10 * (model.x[i] - 0.4) ** 2 for i in range(2, n+m))).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[1, 0] = 0
+    model.Df[0, 1] = 0
+    model.Df[1, 1] = 1
+    for j in range(2,n+m):
+        model.Df[0, j] = 0
+        model.Df[1, j] = 20 * model.x[j] - 8
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = sum(model.x[i] ** 2 for i in range(n+m)) - 4 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    for j in range(n+m):
+        model.Dg[0, j] = 2 * model.x[j]
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-2,-2]])
+    model_parameters.U = np.array([[2,2+(n+m-2)*1.6]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T4.py b/python/problems/T4.py
new file mode 100644
index 0000000..b3e7d3d
--- /dev/null
+++ b/python/problems/T4.py
@@ -0,0 +1,83 @@
+"""
+T4 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+
+It is taken from:
+Marianna De Santis, Gabriele Eichfelder, Julia Niebling, Stefan
+Rocktäschel. Solving Multiobjective Mixed Integer Convex Optimization
+Problems, SIAM Journal on Optimization (2020)
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(params):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = params[0] # Continuous variables
+    m = params[1] # Integer variables
+    p = 2         # Dimension criterion space
+    q = 1         # Number of constraints
+    assert n % 2 == 0, 'Number of continuous variables has to be even.'
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-2,2)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = sum(model.x[i] for i in range(n // 2)) + sum(model.x[i] for i in range(n, n + m))).deactivate()
+    model.f.add(expr = sum(model.x[i] for i in range(n // 2, n)) - sum(model.x[i] for i in range(n, n + m))).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    for j in range(n//2):
+        model.Df[0, j] = 1
+        model.Df[1, j] = 0
+    for j in range(n//2, n):
+        model.Df[0, j] = 0
+        model.Df[1, j] = 1
+    for j in range(n,n+m):
+        model.Df[0, j] = 1
+        model.Df[1, j] = -1
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = sum(model.x[i] ** 2 for i in range(n)) - 1 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    for j in range(n):
+        model.Dg[0, j] = 2 * model.x[j]
+    for j in range(n, n+m):
+        model.Dg[0, j] = 0
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-2*m-n/2,-2*m-n/2]])
+    model_parameters.U = np.array([[2*m+n/2,2*m+n/2]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T5.py b/python/problems/T5.py
new file mode 100644
index 0000000..93bcdc6
--- /dev/null
+++ b/python/problems/T5.py
@@ -0,0 +1,88 @@
+"""
+T5 is a tri-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+
+It is taken from:
+Marianna De Santis, Gabriele Eichfelder, Julia Niebling, Stefan
+Rocktäschel. Solving Multiobjective Mixed Integer Convex Optimization
+Problems, SIAM Journal on Optimization (2020)
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 3           # Continuous variables
+    m = 1           # Integer variables
+    p = 3           # Dimension criterion space
+    q = 1           # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-2,2)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[3]).deactivate()
+    model.f.add(expr = model.x[1] - model.x[3]).deactivate()
+    model.f.add(expr = model.x[2] + model.x[3] ** 2).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 0
+    model.Df[0, 2] = 0
+    model.Df[0, 3] = 1
+
+    model.Df[1, 0] = 0
+    model.Df[1, 1] = 1
+    model.Df[1, 2] = 0
+    model.Df[1, 3] = -1
+
+    model.Df[2, 0] = 0
+    model.Df[2, 1] = 0
+    model.Df[2, 2] = 1
+    model.Df[2, 3] = 2 * model.x[3]
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = sum(model.x[i] ** 2 for i in range(n)) - 1 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    for j in range(n):
+        model.Dg[0, j] = 2 * model.x[j]
+    for j in range(n, n+m):
+        model.Dg[0, j] = 0
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-3,-3,-1]])
+    model_parameters.U = np.array([[3,3,5]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T6.py b/python/problems/T6.py
new file mode 100644
index 0000000..9c00361
--- /dev/null
+++ b/python/problems/T6.py
@@ -0,0 +1,78 @@
+"""
+T6 is a bi-objective mixed-integer convex test instance with a
+non-quadratic objective function and quadratic constraint functions.
+
+It is taken from:
+Marianna De Santis, Gabriele Eichfelder, Julia Niebling, Stefan
+Rocktäschel. Solving Multiobjective Mixed Integer Convex Optimization
+Problems, SIAM Journal on Optimization (2020)
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 2         # Continuous variables
+    m = 1         # Integer variables
+    p = 2         # Dimension criterion space
+    q = 1         # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = False
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-2,2)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-2,2)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[2]).deactivate()
+    model.f.add(expr = model.x[1] + exp(-model.x[2])).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 0
+    model.Df[0, 2] = 1
+    model.Df[1, 0] = 0   
+    model.Df[1, 1] = 1
+    model.Df[1, 2] = exp(-model.x[2])
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = model.x[0] ** 2 + model.x[1] ** 2 - 1 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    model.Dg[0, 0] = 2 * model.x[0]
+    model.Dg[0, 1] = 2 * model.x[1]
+    model.Dg[0, 2] = 0
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-3,-1+np.exp(-2)]])
+    model_parameters.U = np.array([[3,1+np.exp(2)]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T8.py b/python/problems/T8.py
new file mode 100644
index 0000000..08defde
--- /dev/null
+++ b/python/problems/T8.py
@@ -0,0 +1,119 @@
+"""
+T8 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 4         # Continuous variables
+    m = 4         # Integer variables
+    p = 2         # Dimension criterion space
+    q = 4         # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-20,20)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (0,20)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[2] ** 2 + model.x[4] + model.x[6] ** 2).deactivate()
+    model.f.add(expr = model.x[1] ** 2 + model.x[3] + model.x[5] ** 2 + model.x[7]).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 0
+    model.Df[0, 2] = 2 * model.x[2]
+    model.Df[0, 3] = 0
+    model.Df[0, 4] = 1
+    model.Df[0, 5] = 0
+    model.Df[0, 6] = 2 * model.x[6]
+    model.Df[0, 7] = 0
+
+    model.Df[1, 0] = 0
+    model.Df[1, 1] = 2 * model.x[1]
+    model.Df[1, 2] = 0
+    model.Df[1, 3] = 1
+    model.Df[1, 4] = 0
+    model.Df[1, 5] = 2 * model.x[5]
+    model.Df[1, 6] = 0
+    model.Df[1, 7] = 1
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = model.x[0] ** 2 + model.x[1] ** 2 - 1 <= 0)
+    model.g.add(expr = (model.x[2] - 2) ** 2 + (model.x[3] - 5) ** 2 - 10 <= 0)
+    model.g.add(expr = model.x[4] + model.x[5] + model.x[6] + model.x[7] - 5 <= 0)
+    model.g.add(expr = (model.x[4] + model.x[7]) ** 2 - 4 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    model.Dg[0, 0] = 2 * model.x[0]
+    model.Dg[0, 1] = 2 * model.x[1]
+    model.Dg[0, 2] = 0
+    model.Dg[0, 3] = 0
+    model.Dg[0, 4] = 0
+    model.Dg[0, 5] = 0
+    model.Dg[0, 6] = 0
+    model.Dg[0, 7] = 0
+
+    model.Dg[1, 0] = 0
+    model.Dg[1, 1] = 0
+    model.Dg[1, 2] = 2 * (model.x[2] - 2)
+    model.Dg[1, 3] = 2 * (model.x[3] - 5)
+    model.Dg[1, 4] = 0
+    model.Dg[1, 5] = 0
+    model.Dg[1, 6] = 0
+    model.Dg[1, 7] = 0
+
+    model.Dg[2, 0] = 0
+    model.Dg[2, 1] = 0
+    model.Dg[2, 2] = 0
+    model.Dg[2, 3] = 0
+    model.Dg[2, 4] = 1
+    model.Dg[2, 5] = 1
+    model.Dg[2, 6] = 1
+    model.Dg[2, 7] = 1
+
+    model.Dg[3, 0] = 0
+    model.Dg[3, 1] = 0
+    model.Dg[3, 2] = 0
+    model.Dg[3, 3] = 0
+    model.Dg[3, 4] = 2 * (model.x[4] + model.x[7])
+    model.Dg[3, 5] = 0
+    model.Dg[3, 6] = 0
+    model.Dg[3, 7] = 2 * (model.x[4] + model.x[7])
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.empty((1,p))
+    model_parameters.U = np.empty((1,p))
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/problems/T9.py b/python/problems/T9.py
new file mode 100644
index 0000000..23ff672
--- /dev/null
+++ b/python/problems/T9.py
@@ -0,0 +1,119 @@
+"""
+T9 is a bi-objective mixed-integer convex test instance with
+quadratic objective and constraint functions.
+"""
+
+import numpy as np
+from pyomo.environ import *
+
+class structure():
+    pass
+
+def build_model(_):
+
+    model = ConcreteModel()
+    model_parameters = structure()
+
+    # Dimension of decision and criterion space
+    n = 4         # Continuous variables
+    m = 4         # Integer variables
+    p = 2         # Dimension criterion space
+    q = 4         # Number of constraints
+
+    model_parameters.n = n
+    model_parameters.m = m
+    model_parameters.p = p
+    model_parameters.q = q
+
+    # Problem type
+    model_parameters.is_convex = True
+    model_parameters.is_quadratic = True
+
+    # Variables
+    model.x = Var(range(n+m))
+    for i in range(n):
+        model.x[i].domain = Reals
+        model.x[i].bounds = (-20,20)
+    for j in range(n,n+m):
+        model.x[j].domain = Integers
+        model.x[j].bounds = (-20,20)
+
+    # Objective functions
+    model.f = ObjectiveList()
+    model.f.add(expr = model.x[0] + model.x[2] + model.x[4] + model.x[6]).deactivate()
+    model.f.add(expr = model.x[1] + model.x[3] + model.x[5] + model.x[7]).deactivate()
+
+    # Derivatives of objective functions
+    model.Df = Expression(range(p), range(n+m))
+    model.Df[0, 0] = 1
+    model.Df[0, 1] = 0
+    model.Df[0, 2] = 1
+    model.Df[0, 3] = 0
+    model.Df[0, 4] = 1
+    model.Df[0, 5] = 0
+    model.Df[0, 6] = 1
+    model.Df[0, 7] = 0
+
+    model.Df[1, 0] = 0
+    model.Df[1, 1] = 1
+    model.Df[1, 2] = 0
+    model.Df[1, 3] = 1
+    model.Df[1, 4] = 0
+    model.Df[1, 5] = 1
+    model.Df[1, 6] = 0
+    model.Df[1, 7] = 1
+
+    # Constraint functions
+    model.g = ConstraintList()
+    model.g.add(expr = model.x[0] ** 2 + model.x[1] ** 2 - 1 <= 0)
+    model.g.add(expr = model.x[2] ** 2 + model.x[3] ** 2 - 1 <= 0)
+    model.g.add(expr = (model.x[4] - 2) ** 2 + (model.x[5] - 5) ** 2 - 10 <= 0)
+    model.g.add(expr = (model.x[6] - 3) ** 2 + (model.x[7] - 8) ** 2 - 10 <= 0)
+
+    # Derivatives of constraint functions
+    model.Dg = Expression(range(q), range(n+m))
+    model.Dg[0, 0] = 2 * model.x[0]
+    model.Dg[0, 1] = 2 * model.x[1]
+    model.Dg[0, 2] = 0
+    model.Dg[0, 3] = 0
+    model.Dg[0, 4] = 0
+    model.Dg[0, 5] = 0
+    model.Dg[0, 6] = 0
+    model.Dg[0, 7] = 0
+
+    model.Dg[1, 0] = 0
+    model.Dg[1, 1] = 0
+    model.Dg[1, 2] = 2 * model.x[2]
+    model.Dg[1, 3] = 2 * model.x[3]
+    model.Dg[1, 4] = 0
+    model.Dg[1, 5] = 0
+    model.Dg[1, 6] = 0
+    model.Dg[1, 7] = 0
+
+    model.Dg[2, 0] = 0
+    model.Dg[2, 1] = 0
+    model.Dg[2, 2] = 0
+    model.Dg[2, 3] = 0
+    model.Dg[2, 4] = 2 * (model.x[4] - 2)
+    model.Dg[2, 5] = 2 * (model.x[5] - 5)
+    model.Dg[2, 6] = 0
+    model.Dg[2, 7] = 0
+
+    model.Dg[3, 0] = 0
+    model.Dg[3, 1] = 0
+    model.Dg[3, 2] = 0
+    model.Dg[3, 3] = 0
+    model.Dg[3, 4] = 0
+    model.Dg[3, 5] = 0
+    model.Dg[3, 6] = 2 * (model.x[6] - 3)
+    model.Dg[3, 7] = 2 * (model.x[7] - 8)
+
+    # Solver
+    model_parameters.solver = 'gurobi'
+    model_parameters.patch_solver = 'ipopt'
+
+    # Initial enclosure
+    model_parameters.L = np.array([[-3,5]])
+    model_parameters.U = np.array([[13,22]])
+
+    return model, model_parameters
\ No newline at end of file
diff --git a/python/requirements.txt b/python/requirements.txt
new file mode 100644
index 0000000..0b7cf8c
--- /dev/null
+++ b/python/requirements.txt
@@ -0,0 +1,6 @@
+gurobi
+ipopt
+notebook
+numpy
+matplotlib
+pyomo
\ No newline at end of file
diff --git a/python/solver/HyPaD.py b/python/solver/HyPaD.py
new file mode 100644
index 0000000..29d54d9
--- /dev/null
+++ b/python/solver/HyPaD.py
@@ -0,0 +1,282 @@
+import numpy as np
+import time
+import warnings
+from solver.update_bounds import *
+from solver.subsolvers import *
+from solver.snia_functions import *
+
+class IdsEntry:
+    def __init__(self, x_int, L, active_state, nd_candidates, eff_candidates):
+        self.x_int = x_int
+        self.L = L
+        self.S = active_state
+        self.N = nd_candidates
+        self.E = eff_candidates
+
+def f_eval(model, x):
+    """
+    Evaluates the objective functions of the model at a certain point
+
+    Parameters
+    ----------
+    model : pyomo model
+        a model respresenting the current optimization problem
+    x : numpy array
+        point where to evaluate the objective functions
+
+    Returns
+    -------
+    value_obj : numpy array
+        an array consisting of the objective function values
+
+    """
+
+    dim = len(x)
+    #vars = list(model.component_objects(Var))
+    for i in range(dim):
+        #vars[i].value = x[i]
+        model.x[i].value = x[i]
+    #value_obj = np.array([obj.expr() for obj in model.component_objects(Objective)])
+    value_obj = np.array([model.f[i+1].expr() for i in range(len(model.f))])
+
+    return value_obj
+
+def compute_nds(list):
+    """
+    Computes the nondominated points of a given list
+
+    Parameters
+    ----------
+    list : numpy array
+        array of all points
+
+    Returns
+    -------
+    indexlist : numpy array
+        array of boolean values indicating which of the points in the list are nondominated (True) or dominated (False)
+
+    """
+
+    size_list = list.shape[0]
+    indexlist = np.ones(size_list, dtype=bool)
+
+    for i in range(size_list):
+        current_point = list[i,:]
+        indexlist[i] = False
+        if not np.any(np.all(current_point >= list[indexlist,:], axis = 1)):
+            indexlist[i] = True
+    
+    return indexlist
+
+def HyPaD(build_model, params, L, U, EPSILON, OFFSET):
+    """
+    Main procedure for the Hybrid Patch Decomposition algorithm (HyPaD)
+
+    Parameters
+    ----------
+    build_model : function
+        function that returns the pyomo model of the optimization problem and the model paramaters
+    params : array
+        parameters for the model of the optimization problem
+    L : numpy array
+        initial array of lower bounds for the enclosure (with each lower bound being a numpy array itself)
+    U : numpy array
+        initial array of upper bounds for the enclosure (with each upper bound being a numpy array itself)
+    EPSILON : float
+        tolerance for the width of the enclosure
+    OFFSET : float
+        tolerance for comparisons and the evaluation of inequalities
+
+    Returns
+    -------
+    L : numpy array
+        array of lower bounds for the enclosure (with each lower bound being a numpy array itself)
+    U : numpy array
+        array of upper bounds for the enclosure (with each upper bound being a numpy array itself)
+    N : numpy array
+        array of all potentially nondominated points computed by the algorithm
+    ids : array
+        integer data structure containing the data computed on the patch level
+    it : int
+        number of iterations
+    exitflag : int
+        value that indicates why the algorithm terminated
+
+    """
+    
+    # Load model parameters
+    model,model_parameters = build_model(params)
+    n = model_parameters.n
+    m = model_parameters.m
+    p = model_parameters.p
+
+    # Initialization phase    
+    N = np.empty((0,p))
+    E = np.empty((0,n+m))
+    X = np.empty((0,n+m))
+    ids = []
+    it = 0
+    done = False
+    exitflag = -3
+    TIME_LIMIT = 3600
+    ALL_INTEGERS_VISITED = False
+    boxesLB,boxesUB,int_num,current_int = snia_boxes_init(m, np.array([model.x[i].lb for i in range(n,n+m)]), np.array([model.x[i].ub for i in range(n,n+m)]), 4)
+
+    # Test a first solution
+    l = np.min(L, axis = 0)
+    r = np.max(U, axis = 0) - l
+    gsup_x, _, _ = GSUP(build_model, params, l, r)
+    gsup_x_int = np.round(gsup_x[n:n+m])
+    x, err = FeasibilityCheckPatch(build_model, params, gsup_x_int)
+    if err < OFFSET:
+        eff_candidates = WSNLP(build_model, params,gsup_x_int, np.eye(p))
+        nd_candidates = np.array([f_eval(model,eff_candidates[i,:]) for i in range(eff_candidates.shape[0])])
+        z = np.min(nd_candidates, axis = 0)
+        N = np.row_stack([N, nd_candidates])
+        E = np.row_stack([E, eff_candidates])
+        for k in range(nd_candidates.shape[0]):
+            U = updateLUB3(U,nd_candidates[k,:])
+        ids.append(IdsEntry(gsup_x_int,np.array([z-OFFSET]),True,nd_candidates,eff_candidates))
+    else:
+        X = np.row_stack([X, x])
+    
+    # Main Loop
+    start_time = time.time()
+    while (time.time() - start_time < TIME_LIMIT) and not done:
+        it += 1
+        size_L = L.shape[0]
+        finished_bounds = 0
+        rsup_error = True
+        for l in L:
+            if (time.time() - start_time > TIME_LIMIT):
+                break  # Recognize time limit early
+            U_temp = np.all(l < (U - EPSILON + OFFSET), axis=1)    
+            if np.any(U_temp):
+                directions = U[U_temp,:] - l
+                u_index = np.argmax(np.min(directions, axis=1))
+                r = directions[u_index,:]
+                u = l + r
+                rsup_x, rsup_eta, rsup_flag = RSUP(build_model, params, np.vstack((X, E)), l, r) 
+                if rsup_flag < 0:
+                    warnings.warn("RSUP terminated with non-optimal solution", category=UserWarning)
+                else:
+                    rsup_x_int = np.round(rsup_x[n:n+m])
+                    rsup_error = False
+                    if np.all(rsup_eta > u + OFFSET):
+                        # Try to improve the lower bound set using continuous relaxation
+                        warnings.warn("RSUP returned eta > u", category=UserWarning)
+                        gsup_x, gsup_t, gsup_flag = GSUP(build_model, params, l, r)
+                        if gsup_flag < 0:
+                            warnings.warn("GSUP terminated with non-optimal solution", category=UserWarning)
+                        else:
+                            L = updateLLB3(L, l + gsup_t * r)
+                    else:
+                        L = updateLLB3(L, rsup_eta)
+                    if ~ALL_INTEGERS_VISITED:
+                        x, err = FeasibilityCheckPatch(build_model, params, rsup_x_int)
+                        if err < OFFSET:
+                            if any(ids):
+                                ids_index = next((i for i, entry in enumerate(ids) if np.all(entry.x_int == rsup_x_int)), None)
+                            else:
+                                ids_index = None
+                            if ids_index == None:
+                                # InitIDS: Initialize entry of ids
+                                eff_candidates = WSNLP(build_model, params,rsup_x_int, np.eye(p))
+                                nd_candidates = np.array([f_eval(model,eff_candidates[i,:]) for i in range(eff_candidates.shape[0])])
+                                z = np.min(nd_candidates, axis = 0)
+                                N = np.row_stack([N, nd_candidates])
+                                E = np.row_stack([E, eff_candidates])
+                                for k in range(nd_candidates.shape[0]):
+                                    U = updateLUB3(U,nd_candidates[k,:])
+                                ids.append(IdsEntry(rsup_x_int,np.array([z-OFFSET]),True,nd_candidates,eff_candidates))
+                            else:
+                                if not ids[ids_index].S:
+                                    # Select another active patch since the current one is inactive
+                                    ids_index = next((i for i, entry in enumerate(ids) if entry.S), None)
+                                if not ids_index == None:
+                                    # UpdateIDS: Update entry of IDS
+                                    patch_L = ids[ids_index].L
+                                    patch_done = True
+                                    for patch_l in patch_L:
+                                        if (time.time() - start_time > TIME_LIMIT):
+                                            break  # Recognize time limit early
+                                        patch_U_temp = np.all(patch_l < (U - EPSILON + OFFSET), axis=1)
+                                        
+                                        if np.any(patch_U_temp):
+                                            patch_done = False
+                                            patch_directions = U[patch_U_temp,:] - patch_l
+                                            patch_u_index = np.argmax(np.min(patch_directions, axis=1))
+                                            patch_r = patch_directions[patch_u_index,:]
+                                            
+                                            sup_x, sup_t, sup_flag = SUP(build_model, params, ids[ids_index].x_int, patch_l, patch_r)
+                                            
+                                            if sup_flag < 0:
+                                                warnings.warn("SUP terminated with non-optimal solution", category=UserWarning)
+                                            else:
+                                                y_lub = f_eval(model,sup_x)
+                                                y_llb = patch_l + sup_t * patch_r
+                                                U = updateLUB3(U, y_lub)
+                                                ids[ids_index].L = updateLLB3(ids[ids_index].L, y_llb)
+                                                ids[ids_index].N = np.row_stack([ids[ids_index].N, y_lub])
+                                                ids[ids_index].E = np.row_stack([ids[ids_index].E, sup_x])
+                                                N = np.row_stack([N, y_lub])
+                                                E = np.row_stack([E, sup_x])
+                                    if patch_done:
+                                        ids[ids_index].S = False
+                                else:
+                                    # SNIA: Compute new integer assignment
+                                    #TODO: Implement more SNIA modes
+                                    snia_x,snia_flag,ALL_INTEGERS_VISITED,current_int = snia_boxes_fixed(build_model, params, ids, X, boxesLB, boxesUB, int_num, current_int, l, r, OFFSET)
+                                    if snia_flag > 0:
+                                        # InitIDS: Initialize entry of ids
+                                        snia_x_int = snia_x[n:n+m]
+                                        eff_candidates = WSNLP(build_model, params,snia_x_int, np.eye(p))
+                                        nd_candidates = np.array([f_eval(model,eff_candidates[i,:]) for i in range(eff_candidates.shape[0])])
+                                        z = np.min(nd_candidates, axis = 0)
+                                        N = np.row_stack([N, nd_candidates])
+                                        E = np.row_stack([E, eff_candidates])
+                                        for k in range(nd_candidates.shape[0]):
+                                            U = updateLUB3(U,nd_candidates[k,:])
+                                        ids.append(IdsEntry(snia_x_int,np.array([z-OFFSET]),True,nd_candidates,eff_candidates))
+                                        #print('INFO: SNIA computed a new feasible integer assignment.')
+                                    else:
+                                        if ALL_INTEGERS_VISITED:
+                                            # Terminate HyPaD and compute global lower bound set based on patch level lower bounds
+                                            print('END:  All integer assigments visited! Switching solver mode.')
+                                            L = np.vstack([entry.L for entry in ids])
+                                            L_indexlist = compute_nds(L)
+                                            L = L[L_indexlist,:]
+                                            exitflag = 0
+                                            done = True
+                                            return L, U, N, ids, it, exitflag
+                                        else:
+                                            if snia_flag < 0:
+                                                print('INFO: SNIA computed an invalid integer assignment.')
+                                                pass
+                                            else:
+                                                X = np.row_stack([X, snia_x])
+                                                print('INFO: SNIA computed a new infeasible integer assignment.')
+                        else:
+                            X = np.row_stack([X, x])
+                indexlist = [i for i, entry in enumerate(ids) if entry.S]
+            else:
+                finished_bounds += 1
+        if finished_bounds > size_L-1:
+            print("END:  Finished! Width is <= EPSILON.")
+            exitflag = 1
+            done = True
+            break
+        if rsup_error:
+            print("ERR:  Algorithm terminated due to RSUP issues.")
+            exitflag = -2
+            done = True
+            break
+        if not np.any(indexlist):
+            #print('INFO: Iteration ended with all patches inactive');
+            pass
+    if time.time() - start_time > TIME_LIMIT:
+        exitflag = -1
+    elif done:
+        exitflag = 1
+    
+    return L, U, N, ids, it, exitflag
\ No newline at end of file
diff --git a/python/solver/plot_functions.py b/python/solver/plot_functions.py
new file mode 100644
index 0000000..526db97
--- /dev/null
+++ b/python/solver/plot_functions.py
@@ -0,0 +1,106 @@
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+from mpl_toolkits.mplot3d import Axes3D
+import numpy as np
+
+def plot_boxes(L, U, p):
+    """
+    Plots the box coverage given a set of lower and upper bounds for the two- and three-dimensional setting
+
+    Parameters
+    ----------
+    L : numpy array
+        an array of lower bounds (with each lower bound being a numpy array itself)
+    U : numpy array
+        an array of upper bounds (with each upper bound being a numpy array itself)
+    p : int
+        dimension of the boxes
+
+    """
+
+    if p < 3:
+        fig = plt.figure()
+        ax = fig.add_subplot(111)
+        ax.set_xlabel('$f_1$', fontsize=15)
+        ax.set_ylabel('$f_2$', fontsize=15)
+        ax.grid(True)
+        
+        vert_base = np.array([[0, 0, 1, 1], [0, 1, 1, 0]])
+        faces = [1, 2, 4, 3, 1]
+        
+        for l in L:
+            U_boxes = np.all(l < U + 1e-6, axis=1)
+            if np.any(U_boxes):
+                directions = U[U_boxes,:] - l
+                for r in directions:
+                    ax.fill(l[0] + vert_base[0,:] * r[0], l[1] + vert_base[1,:] * r[1],
+                            facecolor=[0, 0.455, 0.478], alpha=0.1, edgecolor=[0, 0.455, 0.478])
+    
+    elif p < 4:
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        ax.set_xlabel('$f_1$', fontsize=15)
+        ax.set_ylabel('$f_2$', fontsize=15)
+        ax.set_zlabel('$f_3$', fontsize=15)
+        ax.grid(True)
+        ax.view_init(elev=15, azim=-130)
+        
+        vert_base = np.array([[0, 0, 0, 0, 1, 1, 1, 1],
+                              [0, 0, 1, 1, 0, 0, 1, 1],
+                              [0, 1, 0, 1, 0, 1, 0, 1]])
+        faces = np.array([[1, 2, 4, 3], [3, 4, 8, 7], [7, 8, 6, 5],
+                          [5, 7, 3, 1], [1, 5, 6, 2], [2, 4, 8, 6]])
+        
+        for l in L:
+            U_boxes = np.all(l < U + 1e-3, axis=1)
+            if np.any(U_boxes):
+                directions = U[U_boxes,:] - l + 1e-6
+                for r in directions:
+                    verts = np.array([l[0] + vert_base[0,:] * r[0],
+                                      l[1] + vert_base[1,:] * r[1],
+                                      l[2] + vert_base[2,:] * r[2]])
+                    ax.add_collection3d(Poly3DCollection(np.array([verts[:,face-1].T for face in faces]),
+                                                        facecolors=[0, 0.455, 0.478], alpha=0.1, edgecolors=[0, 0.455, 0.478]))
+        
+        plot_lb = np.min(L,axis=0)-1
+        plot_ub = np.max(U,axis=0)+1
+        ax.set_xlim(plot_lb[0], plot_ub[0])
+        ax.set_ylim(plot_lb[1], plot_ub[1])
+        ax.set_zlim(plot_lb[2], plot_ub[2])
+    
+    plt.show()
+
+def plot_bounds(L, U, p):
+    """
+    Plots a set of lower and a set of upper bounds for the two- and three-dimensional setting
+
+    Parameters
+    ----------
+    L : numpy array
+        an array of lower bounds (with each lower bound being a numpy array itself)
+    U : numpy array
+        an array of upper bounds (with each upper bound being a numpy array itself)
+    p : int
+        dimension of the bounds
+
+    """
+    
+    if p < 3:
+        fig = plt.figure()
+        ax = fig.add_subplot(111)
+        ax.plot(L[:,0], L[:,1], linestyle='none', marker='.', color='blue')
+        ax.plot(U[:,0], U[:,1], linestyle='none', marker='.', color=[1, 0.4745, 0])
+        ax.grid(True)
+        ax.set_xlabel('$f_1$')
+        ax.set_ylabel('$f_2$')
+    elif p < 4:
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        ax.plot(L[:,0], L[:,1], L[:,2], linestyle='none', marker='.', color='blue')
+        ax.plot(U[:,0], U[:,1], U[:,2], linestyle='none', marker='.', color=[1, 0.4745, 0])
+        ax.grid(True)
+        ax.set_xlabel('$f_1$')
+        ax.set_ylabel('$f_2$')
+        ax.set_zlabel('$f_3$')
+        ax.view_init(elev=15, azim=-130)
+    plt.show()
\ No newline at end of file
diff --git a/python/solver/snia_functions.py b/python/solver/snia_functions.py
new file mode 100644
index 0000000..b234f1c
--- /dev/null
+++ b/python/solver/snia_functions.py
@@ -0,0 +1,154 @@
+import numpy as np
+from solver.subsolvers import *
+
+def snia_boxes_init(m, int_lb, int_ub, num_splits):
+    """
+    Initializes the boxes to search for a new integer assignment for the fixed boxes approach
+
+    Parameters
+    ----------
+    m : int
+        number of integer variables
+    int_lb : numpy array
+        lower bounds for the integer variables
+    int_ub : numpy array
+        upper bounds for the integer variables
+    num_splits : int
+        number of splits
+
+    Returns
+    -------
+    boxesLB : numpy array
+        lower bounds of all boxes (with each lower bound being a numpy array itself)
+    boxesUB : numpy array
+        upper bounds of all boxes (with each upper bound being a numpy array itself)
+    int_num : int
+        number of possible integer assignments
+    current_int : numpy array
+        array that indicates the position of the last integer assignment that has been computed for each box (which is 0 after the initialization)
+
+    """
+
+    # Extracting integer bounds
+    width = int_ub - int_lb
+    int_num = np.prod(width + 1)
+    boxesLB = np.array([int_lb])
+    boxesUB = np.array([int_ub])
+
+    for i in range(num_splits):
+        max_width, k = np.max(width), np.argmax(width)
+        if max_width < 0.75:
+            break
+        
+        offset = np.zeros(m)
+        offset[k] = np.ceil(max_width / 2 - 0.25)
+        
+        boxesLB = np.row_stack((boxesLB, boxesLB + offset))
+        boxesUB = np.row_stack((boxesUB - offset, boxesUB))
+        
+        width[k] = np.floor(max_width / 2 + 0.25)
+
+    current_int = np.zeros(boxesLB.shape[0], dtype=int)
+
+    return boxesLB, boxesUB, int_num, current_int
+
+
+def snia_boxes_fixed(build_model, params, ids, X, boxesLB, boxesUB, int_num, current_int, a, r, OFFSET ,GUESS=False):
+    """
+    Performs the search for a new integer assignment based on the fixed boxes approach
+
+    Parameters
+    ----------
+    build_model : function
+        function that returns the pyomo model of the optimization problem and the model paramaters
+    params : array
+        parameters for the model of the optimization problem
+    ids : numpy array
+        point from where to start the search for an update point
+    X : numpy array
+        direction in which to search for an update point
+    boxesLB : numpy array
+        lower bounds of all boxes (with each lower bound being a numpy array itself)
+    boxesUB : numpy array
+        upper bounds of all boxes (with each upper bound being a numpy array itself)
+    int_num : int
+        number of possible integer assignments
+    current_int : numpy array
+        array that indicates the position of the last integer assignment that has been computed for each box
+    a : numpy array
+        parameter for GSUP when using the buest guess based on continuous relaxation
+    r : numpy array
+        parameter for GSUP when using the buest guess based on continuous relaxation
+    OFFSET : float
+        tolerance for comparisons and the evaluation of inequalities
+    GUESS : boolean
+        indicates whether to first try to compute a new integer assignment based on continuous relaxation
+
+    Returns
+    -------
+    sol_x : numpy array
+        part of an optimal solution of SUP corresponding to the x variable of the original optimization problem
+    flag : int
+        value that indicates why the procedure terminated
+    allSolved : boolean
+        indicates whether all integer assignments have been visited
+    current_int : numpy array
+        updated array that indicates the position of the last integer assignment that has been computed for each box
+
+    """
+
+    # Searches a new integer assignment within predefined boxes
+
+    # Load the model parameters
+    model, model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+    q = model_param.q
+
+    # Initialization
+    intArray = np.row_stack((np.row_stack([entry.x_int for entry in ids]), X[:, n:n+m]))
+    allSolved = False
+    flag = -1
+
+    # Check if all integers are already solved
+    if not intArray.shape[0] < int_num:
+        allSolved = True
+        sol_x = np.array([model.x[i].lb for i in model.x])
+        flag = 0
+        return sol_x, flag, allSolved, current_int
+
+    # Buest guess based on relaxation
+    if GUESS:
+        gsup_x, _, _ = GSUP(build_model, params, a, r)
+        x_int = np.round(gsup_x[n:n+m])
+        sol_x,err = FeasibilityCheckPatch(build_model, params, x_int)
+        if err < OFFSET:
+            flag = 1
+        else:
+            flag = 0
+        if not np.any(np.all(intArray == x_int, axis = 1)):
+            return sol_x, flag, allSolved, current_int
+
+    # Find next integer assignment
+    boxes_num = boxesLB.shape[0]
+    count = np.zeros(boxes_num)
+
+    for i in range(boxes_num):
+        count[i] = np.sum(np.all((boxesLB[i,:] <= intArray) & (intArray <= boxesUB[i,:]), axis=1))
+
+    index = np.argmin(count)
+    int_lb = boxesLB[index, :]
+    int_ub = boxesUB[index, :]
+    width = int_ub - int_lb
+    width_prod = np.prod(np.triu(np.ones((m, m)),1) * (width + 1) + np.tril(np.ones((m, m))), axis=1)
+
+    for i in range(current_int[index], int(np.prod(width + 1))):
+        x_int = int_lb + np.mod(np.floor(i / width_prod), width + 1)
+        current_int[index] = i
+        if not np.any(np.all(intArray == x_int, axis = 1)):
+            sol_x,err = FeasibilityCheckPatch(build_model, params, x_int)
+            if err < OFFSET:
+                flag = 1
+            else:
+                flag = 0
+            return sol_x, flag, allSolved, current_int
\ No newline at end of file
diff --git a/python/solver/subsolvers.py b/python/solver/subsolvers.py
new file mode 100644
index 0000000..3f31d48
--- /dev/null
+++ b/python/solver/subsolvers.py
@@ -0,0 +1,310 @@
+import numpy as np
+from pyomo.environ import *
+from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr
+
+def init_enclosure(model, model_parameters, OFFSET):
+    """
+    Computes an initial enclosure consisting of a single upper and lower bound based on interval arithmetic.
+
+    Parameters
+    ----------
+    model : pyomo model
+        a model respresenting the current optimization problem
+    model_parameters : structure
+        a structure respresenting the model parameters
+    OFFSET : float
+        offset for the lower and upper bounds to ensure that the area of interest is contained in the interior of the enclosure
+
+    Returns
+    -------
+    L : numpy array
+        set of lower bounds for the initial enclosure
+    U : numpy array
+        set of upper bounds for the initial enclosure
+
+    """
+
+    p = model_parameters.p
+    L = model_parameters.L
+    U = model_parameters.U
+    if not (np.any(L) and np.any(U)):
+        for i in range(p):
+            L[0,i], U[0,i] = compute_bounds_on_expr(model.f[i+1].expr)
+    return L - OFFSET, U + OFFSET
+
+def WSNLP(build_model, params, x_int, weights):
+    # Load the original model
+    model,model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+
+    # Fix the integer variables
+    for i in range(m):
+        model.x[n+i].fix(x_int[i])
+
+    # Initialization of return values and number of iterations for weighted sum approach to solve the subproblem
+    iterations = weights.shape[0]
+    sol_x = np.empty((iterations,n+m))
+
+    # Solve the Pyomo model
+    solver = SolverFactory(model_param.patch_solver)
+    for i in range(iterations):
+        model.obj = Objective(expr=sum(weights[i,j] * model.f[j+1].expr for j in range(len(model.f))))
+        results = solver.solve(model, tee=False)
+        sol_x[i,:] = np.array([value(model.x[i]) for i in range(len(model.x))])
+        model.del_component(model.obj)
+
+    return sol_x
+
+def FeasibilityCheckPatch(build_model, params, x_int):
+    # Load the original model
+    model,model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+    q = model_param.q
+
+    # Fix the integer variables
+    for i in range(m):
+        model.x[n+i].fix(x_int[i])
+
+    # Add the additional variable alpha
+    model.add_component('alpha', Var(within=Reals))
+    alpha = model.component('alpha')
+
+    # Modify the constraints
+    for j in range(q):
+        con = model.g[j+1]
+        con.set_value(con.body - alpha <= con.upper)
+        
+    # Objective function
+    model.del_component(model.f)
+    model.obj = Objective(expr=alpha)
+
+    # Solve the Pyomo model
+    solver = SolverFactory(model_param.patch_solver)
+    results = solver.solve(model, tee=False)
+
+    # Extract solution
+    err = value(alpha)
+    sol_x = np.array([value(model.x[i]) for i in range(len(model.x))])
+    if not results.solver.termination_condition == TerminationCondition.optimal:
+        exitflag = -1
+    else:
+        exitflag = 0
+
+    return sol_x, err
+
+def SUP(build_model, params, x_int, a, r):
+    """
+    Performs the search for an update point on the patch level
+
+    Parameters
+    ----------
+    build_model : function
+        function that returns the pyomo model of the optimization problem and the model paramaters
+    params : array
+        parameters for the model of the optimization problem
+    x_int : numpy array
+        integer assignment (that describes the current patch problem)
+    a : numpy array
+        point from where to start the search for an update point
+    r : numpy array
+        direction in which to search for an update point
+
+    Returns
+    -------
+    sol_x : numpy array
+        part of an optimal solution of SUP corresponding to the x variable of the original optimization problem
+    sol_t : float
+        optimal value of the optimization problem SUP
+    exitflag : int
+        value that indicates whether the problem SUP was solved to optimality
+
+    """
+
+    # Load the original model
+    model,model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+    p = model_param.p
+
+    # Fix the integer variables
+    for i in range(m):
+        model.x[n+i].fix(x_int[i])
+
+    # Add the additional variable t
+    model.add_component('t', Var(within=Reals))
+    t = model.component('t')
+
+    # Add the SUP constraints
+    for i in range(p):
+        obj = model.f[i+1]
+        model.add_component('sup_con_'+str(i+1),Constraint(expr=obj.expr - t * r[i] <= a[i]))
+
+    # Objective function
+    model.del_component(model.f)
+    model.obj = Objective(expr=t)
+
+    # TODO: Add warm start initialization for x variables
+
+    # Solve the Pyomo model
+    solver = SolverFactory(model_param.patch_solver)
+    results = solver.solve(model, tee=False)
+
+    # Extract solution
+    sol_x = np.array([value(model.x[i]) for i in range(len(model.x))])  
+    sol_t = value(t)
+
+    if not results.solver.termination_condition == TerminationCondition.optimal:
+        exitflag = -1
+    else:
+        exitflag = 0
+
+    return sol_x, sol_t, exitflag
+
+def RSUP(build_model, params, X, a, r):
+    """
+    Performs the search for an update point on the global level based on the relaxed problem R(X)
+
+    Parameters
+    ----------
+    build_model : function
+        function that returns the pyomo model of the optimization problem and the model paramaters
+    params : array
+        parameters for the model of the optimization problem
+    X : numpy array
+        array of linearization points (with each linearization point being a numpy array itself)
+    a : numpy array
+        point from where to start the search for an update point
+    r : numpy array
+        direction in which to search for an update point
+
+    Returns
+    -------
+    sol_x : numpy array
+        part of an optimal solution of RSUP corresponding to the x variable of the relaxed problem R(X)
+    sol_eta : float
+        part of an optimal solution of RSUP corresponding to the eta variable of the relaxed problem R(X)
+    exitflag : int
+        value that indicates whether the problem RSUP was solved to optimality
+
+    """
+    # Load the original model
+    model,model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+    p = model_param.p
+    q = model_param.q
+
+    # Add the additional variable t
+    model.add_component('t', Var(within=Reals))
+    t = model.component('t')
+
+    # Add eta to shift objectives to constraints
+    model.add_component('eta', Var(range(p),within=Reals))
+    eta = model.component('eta')
+
+    # Add the linearized objective and constraint functions
+    for k in range(X.shape[0]):
+        x_lin = X[k,:]
+        for i in range(n+m):
+            model.x[i].value = x_lin[i]
+        for i in range(p):
+            obj = model.f[i+1]
+            model.add_component('f_lin_'+str(k)+'_'+str(i+1),Constraint(expr=obj.expr() + sum(model.Df[i,j].expr() * (model.x[j] - x_lin[j]) for j in range(len(model.x))) - eta[i] <= 0))
+        for i in range(q):
+            con = model.g[i+1]
+            model.add_component('g_lin_'+str(k)+'_'+str(i+1),Constraint(expr=con.body() + sum(model.Dg[i,j].expr() * (model.x[j] - x_lin[j]) for j in range(len(model.x))) <= con.upper))
+
+    # Add the SUP constraints
+    for i in range(p):
+        model.add_component('sup_con_'+str(i+1),Constraint(expr=eta[i] - t * r[i] <= a[i]))
+
+    # Remove old objective and constraint functions
+    model.del_component(model.f)
+    model.del_component(model.g)
+    
+    # Objective function
+    model.obj = Objective(expr=t)
+
+    # Solve the Pyomo model
+    solver = SolverFactory(model_param.solver)
+    results = solver.solve(model, tee=False)
+
+    # Extract solution
+    sol_x = np.array([value(model.x[i]) for i in range(len(model.x))])
+    sol_eta = np.array([value(eta[i]) for i in range(len(eta))])  
+    sol_t = value(t)
+
+    if not results.solver.termination_condition == TerminationCondition.optimal:
+        exitflag = -1
+    else:
+        exitflag = 0
+
+    return sol_x, sol_eta, exitflag
+
+def GSUP(build_model, params, a, r):
+    """
+    Performs the search for an update point on the global level based on the continuous relaxation
+
+    Parameters
+    ----------
+    build_model : function
+        function that returns the pyomo model of the optimization problem and the model paramaters
+    params : array
+        parameters for the model of the optimization problem
+    a : numpy array
+        point from where to start the search for an update point
+    r : numpy array
+        direction in which to search for an update point
+
+    Returns
+    -------
+    sol_x : numpy array
+        part of an optimal solution of GSUP corresponding to the x variable of the original optimization problem
+    sol_t : float
+        optimal value of the optimization problem GSUP
+    exitflag : int
+        value that indicates whether the problem GSUP was solved to optimality
+
+    """
+
+    # Load the original model
+    model,model_param = build_model(params)
+    n = model_param.n
+    m = model_param.m
+    p = model_param.p
+
+    # Change integer variables to reals
+    for j in range(n,n+m):
+        model.x[j].domain = Reals
+
+    # Add the additional variable t
+    model.add_component('t', Var(within=Reals))
+    t = model.component('t')
+
+    # Add the SUP constraints
+    for i in range(p):
+        obj = model.f[i+1]
+        model.add_component('sup_con_'+str(i+1),Constraint(expr=obj.expr - t * r[i] <= a[i]))
+
+    # Objective function
+    model.del_component(model.f)
+    model.obj = Objective(expr=t)
+
+    # TODO: Add warm start initialization for x variables
+
+    # Solve the Pyomo model
+    solver = SolverFactory(model_param.patch_solver)
+    results = solver.solve(model, tee=False)
+
+    # Extract solution
+    sol_x = np.array([value(model.x[i]) for i in range(len(model.x))])  
+    sol_t = value(t)
+
+    if not results.solver.termination_condition == TerminationCondition.optimal:
+        exitflag = -1
+    else:
+        exitflag = 0
+
+    return sol_x, sol_t, exitflag
diff --git a/python/solver/update_bounds.py b/python/solver/update_bounds.py
new file mode 100644
index 0000000..3ec9e9d
--- /dev/null
+++ b/python/solver/update_bounds.py
@@ -0,0 +1,162 @@
+"""
+This file contains the functions to update the sets of local upper and local lower bounds.
+
+Both procedures (especially updateLUB3) implement Algorithm 3 from:
+K. Klamroth et al., On the representation of the search region in multi-objective optimization,
+European Journal of Operational Research (2015)
+"""
+
+import numpy as np
+
+def updateLUB3(ListLUB, z, tolerance=1e-6):
+    """
+    Updates a list of local upper bounds given a new update point
+
+    Parameters
+    ----------
+    ListLUB : numpy array
+        an array of local upper bounds (with each local upper bound being a numpy array itself)
+    z : numpy array
+        update point
+    tolerance : float
+        (optional) tolerance for comparisons when checking inequalitites
+
+    Returns
+    -------
+    ListLUBnew : numpy array
+        updated array of local upper bounds
+
+    """
+
+    # Initialization
+    p = z.shape[0]
+    
+    # A is the set of all search zones {u}-K that contain z
+    A_indexlist = np.all(np.less(z, ListLUB), axis=1)
+    A = ListLUB[A_indexlist,:]
+    sizeA = A.shape[0]
+    
+    B = [None] * p
+    P = [None] * p
+    
+    # Update procedure
+    B[0] = ListLUB[np.all(np.vstack([
+        np.equal(z[0], ListLUB[:,0]).T,
+        (np.less(z[1:p], ListLUB[:,1:p]).T)
+    ]),axis=0),:]
+    P[0] = np.hstack([z[0] * np.ones((sizeA, 1)), A[:,1:p]])
+    for j in range(0,p-1):
+        # B[j] contains all search zones whose boundary contains z with respect
+        # to the j-th component (u_j = z_j)
+        B[j] = ListLUB[np.all(np.vstack([
+            np.less(z[0:j], ListLUB[:, 0:j]).T,
+            np.equal(z[j], ListLUB[:, j]),
+            np.less(z[j+1:p], ListLUB[:, j+1:p]).T
+        ]), axis=0),:]
+        
+        # P[j] contains the projections of z on all local upper bounds in A
+        # along the j-th dimension / component
+        P[j] = np.hstack([A[:, 0:j], z[j] * np.ones((sizeA,1)), A[:,j+1:p]])
+    B[p-1] = ListLUB[np.all(np.vstack([
+        np.less(z[0:p-1], ListLUB[:, 0:p-1]).T,
+        np.equal(z[p-1], ListLUB[:, p-1]),
+    ]), axis=0),:]
+    P[p-1] = np.hstack([A[:, 0:p-1], z[p-1] * np.ones((sizeA,1))])
+    
+    # The following loop filters all redundant points out of P(j) for every
+    # j=1:p where redundant basically means dominated
+    P_indexlist = np.zeros((p, sizeA), dtype=bool)
+    for j in range(p):
+        Pj = P[j]
+        PB = np.vstack([Pj, B[j]])
+        for i in range(sizeA):
+            P_indexlist[j, i] = not np.any(np.all(np.vstack([
+                np.all(np.less_equal(Pj[i, :], PB), axis=1),
+                np.any(np.less(Pj[i, :], PB), axis=1)
+            ]), axis=0))
+        P[j] = Pj[P_indexlist[j, :],:]
+    
+    # P is transformed from a list to a matrix
+    P = np.vstack(P[0:p])
+    
+    # The new list of local upper bounds is computed and made unique
+    ListLUBnew = np.unique(np.vstack([ListLUB[~A_indexlist,:], P]), axis=0)
+    
+    return ListLUBnew
+
+
+def updateLLB3(ListLLB, z, tolerance=1e-6):
+    """
+    Updates a list of local lower bounds given a new update point
+
+    Parameters
+    ----------
+    ListLLB : numpy array
+        an array of local lower bounds (with each local lower bound being a numpy array itself)
+    z : numpy array
+        update point
+    tolerance : float
+        (optional) tolerance for comparisons when checking inequalitites
+
+    Returns
+    -------
+    ListLLBnew : numpy array
+        updated array of local lower bounds
+
+    """
+    
+    # Initialization
+    p = z.shape[0]
+    
+    # A is the set of all search zones {l}+K that contain z
+    A_indexlist = np.all(np.greater(z, ListLLB), axis=1)
+    A = ListLLB[A_indexlist,:]
+    sizeA = A.shape[0]
+    
+    B = [None] * p
+    P = [None] * p
+    
+    # Update procedure
+    B[0] = ListLLB[np.all(np.vstack([
+        np.equal(z[0], ListLLB[:,0]).T,
+        (np.greater(z[1:p], ListLLB[:,1:p]).T)
+    ]),axis=0),:]
+    P[0] = np.hstack([z[0] * np.ones((sizeA, 1)), A[:,1:p]])
+    for j in range(0,p-1):
+        # B[j] contains all search zones whose boundary contains z with respect
+        # to the j-th component (u_j = z_j)
+        B[j] = ListLLB[np.all(np.vstack([
+            np.greater(z[0:j], ListLLB[:, 0:j]).T,
+            np.equal(z[j], ListLLB[:, j]),
+            np.greater(z[j+1:p], ListLLB[:, j+1:p]).T
+        ]), axis=0),:]
+        
+        # P[j] contains the projections of z on all local lower bounds in A
+        # along the j-th dimension / component
+        P[j] = np.hstack([A[:, 0:j], z[j] * np.ones((sizeA,1)), A[:,j+1:p]])
+    B[p-1] = ListLLB[np.all(np.vstack([
+        np.greater(z[0:p-1], ListLLB[:, 0:p-1]).T,
+        np.equal(z[p-1], ListLLB[:, p-1]),
+    ]), axis=0),:]
+    P[p-1] = np.hstack([A[:, 0:p-1], z[p-1] * np.ones((sizeA,1))])
+    
+    # The following loop filters all redundant points out of P(j) for every
+    # j=1:p where redundant basically means dominated
+    P_indexlist = np.zeros((p, sizeA), dtype=bool)
+    for j in range(p):
+        Pj = P[j]
+        PB = np.vstack([Pj, B[j]])
+        for i in range(sizeA):
+            P_indexlist[j, i] = not np.any(np.all(np.vstack([
+                np.all(np.greater_equal(Pj[i, :], PB), axis=1),
+                np.any(np.greater(Pj[i, :], PB), axis=1)
+            ]), axis=0))
+        P[j] = Pj[P_indexlist[j, :],:]
+    
+    # P is transformed from a list to a matrix
+    P = np.vstack(P[0:p])
+    
+    # The new list of local upper bounds is computed and made unique
+    ListLLBnew = np.unique(np.vstack([ListLLB[~A_indexlist,:], P]), axis=0)
+    
+    return ListLLBnew
\ No newline at end of file