Added Python implementation

.gitignore

@@ -1,6 +1,33 @@
-## GitIgnore for HyPaD
-### MATLAB temporary files
-*.aux
-
-### Gurobi problem file 
-*.lp
+## 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

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:

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
+}

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

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