diff --git a/demos/approximate_algorithms.py b/demos/approximate_algorithms.py
new file mode 100644
index 0000000..8f8f605
--- /dev/null
+++ b/demos/approximate_algorithms.py
@@ -0,0 +1,17 @@
+#!python3
+
+"""
+Demonstrates how to use approximation algorithms 
+such as Karmarkar-Karp and Greedy.
+"""
+
+from numberpartitioning import karmarkar_karp, greedy, complete_greedy
+numbers = [4, 6, 7, 5, 8]
+print("Karmarkar-Karp with k=2", karmarkar_karp(numbers, num_parts=2))
+print("Karmarkar-Karp with k=3", karmarkar_karp(numbers, num_parts=3))
+print("Greedy with k=2", greedy(numbers, num_parts=2))
+print("Greedy with k=3", greedy(numbers, num_parts=3))
+
+print("More tests")
+print(greedy([1,2,3,4,5,9], num_parts=2))
+print(greedy([1,2,3,4,5,9], num_parts=3))
diff --git a/demos/complete_algorithms.py b/demos/complete_algorithms.py
new file mode 100644
index 0000000..5b85ffc
--- /dev/null
+++ b/demos/complete_algorithms.py
@@ -0,0 +1,33 @@
+#!python3
+
+"""
+Demonstrates how to use complete algorithms 
+such as Complete Greedy.
+"""
+
+from numberpartitioning import complete_greedy
+
+print("Complete Greedy with k=2, default objective:")
+for p in complete_greedy([4, 6, 7, 5, 8], num_parts=2):
+    print (p)
+
+# Define different objectives (all of them for minimization)
+def largest_sum(partition):
+    return max([sum(part) for part in partition])
+def minus_smallest_sum(partition):
+    return -min([sum(part) for part in partition])
+def largest_diff(partition):
+    sums = [sum(part) for part in partition]
+    return max(sums)-min(sums)
+
+print("An example from Walter (2013), 'Comparing the minimum completion times of two longest-first scheduling-heuristics'.")
+walter_numbers = [46, 39, 27, 26, 16, 13, 10]  
+print("Complete Greedy with k=3, min-diff objective (default):")
+for p in complete_greedy(walter_numbers, num_parts=3):
+    print (p)
+print("Complete Greedy with k=3, min-max objective:")
+for p in complete_greedy(walter_numbers, num_parts=3, objective=largest_sum):
+    print (p)
+print("Complete Greedy with k=3, max-min objective:")
+for p in complete_greedy(walter_numbers, num_parts=3, objective=minus_smallest_sum):
+    print (p)
diff --git a/demos/ilp.py b/demos/ilp.py
new file mode 100644
index 0000000..29a9f60
--- /dev/null
+++ b/demos/ilp.py
@@ -0,0 +1,13 @@
+#!python3
+
+""" 
+A demo program for finding an optimal partition using ILP.
+"""
+
+from numberpartitioning.ilp import *
+
+print("An example from Walter (2013), 'Comparing the minimum completion times of two longest-first scheduling-heuristics'.")
+walter_numbers = [46, 39, 27, 26, 16, 13, 10]  
+print("Min-diff objective:", min_diff_partition(walter_numbers, num_parts=3))
+print("Min-max objective:", min_max_partition(walter_numbers, num_parts=3))
+print("Max-min objective:", max_min_partition(walter_numbers, num_parts=3))
diff --git a/requirements.txt b/requirements.txt
index 1eb213d..aa1c9d1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,5 @@
+cvxpy
+
 # Dev
 build
 
diff --git a/src/numberpartitioning/ilp.py b/src/numberpartitioning/ilp.py
new file mode 100644
index 0000000..6ad16ab
--- /dev/null
+++ b/src/numberpartitioning/ilp.py
@@ -0,0 +1,267 @@
+#!python3
+
+from typing import  List, Optional, Callable
+from numbers import Number
+
+from numberpartitioning.common import Partition, PartitioningResult
+
+import cvxpy
+from numberpartitioning.solve import minimize
+
+
+def max_min_partition(
+    numbers: List[int], num_parts: int = 2, return_indices: bool = False,
+    copies = 1, num_smallest_parts:int = 1,
+) -> PartitioningResult:
+    """
+    Produce a partition that maximizes the smallest sum,
+    by solving an integer linear program (ILP).
+    Credit: Rob Pratt, https://or.stackexchange.com/a/6115/2576
+    Uses CVXPY as an ILP solver.
+
+    Parameters
+    ----------
+    numbers
+        The list of numbers to be partitioned.
+    num_parts
+        The desired number of parts in the partition. Default: 2.
+    return_indices
+        If True, the elements of the parts are the indices of the corresponding entries
+        of numbers; if False (default), the elements are the numbers themselves.
+    copies
+        how many copies are there from each number.
+        Can be either a single integer (the same #copies for all numbers),
+        or a list of integers (a different #copies for each number).
+        Default: 1
+    num_smallest_parts
+        number of smallest parts whose sum should be maximized.
+        Default is 1, which means to just maximize the smallest part-sum.
+        A value of 3, for example, means to maximize the sum of the three smallest part-sums.
+
+    Returns
+    -------
+    A partition represented by a ``PartitioningResult``.
+
+    >>> max_min_partition([10,20,40,0], num_parts=1, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    >>> max_min_partition([10,20,40,0], num_parts=2, return_indices=False)
+    PartitioningResult(partition=[[10, 20, 0], [40]], sizes=[30.0, 40.0])
+
+    >>> result = max_min_partition([10,20,40,0], num_parts=3)
+    >>> int(min(result.sizes))
+    10
+
+    >>> result = max_min_partition([10,20,40,0], num_parts=4)
+    >>> int(min(result.sizes))
+    0
+    
+    >>> result = max_min_partition([10,20,40,0], num_parts=5)
+    >>> int(min(result.sizes))
+    0
+    
+    >>> max_min_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    >>> max_min_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
+
+    >>> max_min_partition([10,20,40,1], num_parts=3, num_smallest_parts=2)
+    PartitioningResult(partition=[[], [10, 20, 1], [40]], sizes=[0.0, 31.0, 40.0])
+    """
+    def minimization_objective(parts_sums):
+        return - sum(parts_sums[0:num_smallest_parts])
+    return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)
+
+
+
+def min_max_partition(
+    numbers: List[int], num_parts: int = 2, return_indices: bool = False,
+    copies = 1, num_largest_parts:int = 1,
+) -> PartitioningResult:
+    """
+    Produce a partition that minimizes the largest sum,
+    by solving an integer linear program (ILP).
+
+    Parameters
+    ----------
+    numbers
+        The list of numbers to be partitioned.
+    num_parts
+        The desired number of parts in the partition. Default: 2.
+    return_indices
+        If True, the elements of the parts are the indices of the corresponding entries
+        of numbers; if False (default), the elements are the numbers themselves.
+    copies
+        how many copies are there from each number.
+        Can be either a single integer (the same #copies for all numbers),
+        or a list of integers (a different #copies for each number).
+        Default: 1
+    num_largest_parts
+        number of largest parts whose sum should be minimized.
+        Default is 1, which means to just minimize the largest part-sum.
+        A value of 3, for example, means to minimize the sum of the three largest part-sums.
+
+    Returns
+    -------
+    A partition represented by a ``PartitioningResult``.
+
+    >>> min_max_partition([10,20,40,0], num_parts=1, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    >>> min_max_partition([10,20,40,0], num_parts=2, return_indices=False)
+    PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])
+
+    >>> result = min_max_partition([10,20,40,0], num_parts=3)
+    >>> int(max(result.sizes))
+    40
+    
+    >>> result = min_max_partition([10,20,40,0], num_parts=5)
+    >>> int(max(result.sizes))
+    40
+    
+    >>> min_max_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    >>> min_max_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
+
+    >>> min_max_partition([10,20,40,1], num_parts=3, num_largest_parts=2)
+    PartitioningResult(partition=[[10, 1], [20], [40]], sizes=[11.0, 20.0, 40.0])
+    """
+    def minimization_objective(parts_sums):
+        return sum(parts_sums[-num_largest_parts:])
+    return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)
+
+
+def min_diff_partition(
+    numbers: List[int], num_parts: int = 2, return_indices: bool = False,
+    copies = 1
+) -> PartitioningResult:
+    """
+    Produce a partition that minimizes the difference between the largest and smallest sum,
+    by solving an integer linear program (ILP).
+
+    Parameters
+    ----------
+    numbers
+        The list of numbers to be partitioned.
+    num_parts
+        The desired number of parts in the partition. Default: 2.
+    return_indices
+        If True, the elements of the parts are the indices of the corresponding entries
+        of numbers; if False (default), the elements are the numbers themselves.
+    copies
+        how many copies are there from each number.
+        Can be either a single integer (the same #copies for all numbers),
+        or a list of integers (a different #copies for each number).
+        Default: 1
+
+    Returns
+    -------
+    A partition represented by a ``PartitioningResult``.
+
+    >>> min_diff_partition([10,20,40,0], num_parts=1, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    >>> min_diff_partition([10,20,40,0], num_parts=2, return_indices=False)
+    PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])
+    
+    >>> min_diff_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    >>> min_diff_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
+    """
+    def minimization_objective(parts_sums):
+        return parts_sums[-1] - parts_sums[0]
+    return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)
+
+
+def find_optimal_partition(
+    numbers: List[int], 
+    minimization_objective: Optional[Callable[[list], float]],
+    num_parts: int = 2, return_indices: bool = False,  copies = 1
+) -> PartitioningResult:
+    """
+    Produce a partition that minimizes the given objective,
+    by solving an integer linear program (ILP).
+    Credit: Rob Pratt, https://or.stackexchange.com/a/6115/2576
+    Uses CVXPY as an ILP solver.
+
+    Parameters
+    ----------
+    numbers
+        The list of numbers to be partitioned.
+    minimization_objective
+        A callable for constructing the objective function to be minimized.
+        Gets as input a list of part-sums, sorted from small to large. 
+        Each of the part-sums is a cvxpy expression.
+        Returns as output an expression that should be minimized.
+        See max_min_partition for usage examples.
+    num_parts
+        The desired number of parts in the partition. Default: 2.
+    return_indices
+        If True, the elements of the parts are the indices of the corresponding entries
+        of numbers; if False (default), the elements are the numbers themselves.
+    copies
+        how many copies are there from each number.
+        Can be either a single integer (the same #copies for all numbers),
+        or a list of integers (a different #copies for each number).
+        Default: 1
+
+    Returns
+    -------
+    A partition representing by a ``PartitioningResult``.
+
+    """
+    parts = range(num_parts)
+    items = range(len(numbers))
+    if isinstance(copies, Number):
+        copies = [copies]*len(numbers)
+
+    counts:dict = {
+        item:
+        [cvxpy.Variable(integer=True) for part in parts]
+        for item in items
+    }	# counts[i][j] determines how many times item i appears in part j.
+    parts_sums = [
+        sum([counts[item][part]*numbers[item] for item in items])
+        for part in parts]
+
+    # Construct the list of constraints:
+    constraints = []
+
+    # The counts must be non-negative:
+    constraints += [counts[item][part]  >= 0 for part in parts for item in items]
+
+    # Each item must be in exactly one part:
+    constraints += [sum([counts[item][part] for part in parts]) == copies[item] for item in items] 	
+
+    # Parts must be in ascending order of their sum (a symmetry-breaker):
+    constraints += [parts_sums[part+1] >= parts_sums[part] for part in range(num_parts-1)]
+
+    objective = minimization_objective(parts_sums)
+    minimize(objective, constraints)
+
+    partition = [
+        sum([int(counts[item][part].value)*[item] 
+            for item in items if counts[item][part].value>=1], [])
+        for part in parts
+    ]
+    sums = [parts_sums[part].value for part in parts]
+    partition:Partition = [[] for _ in parts]
+    for part in parts:
+        for item in items:
+            count = int(counts[item][part].value)
+            if count>=1:
+                partition[part] += count * [item if return_indices else numbers[item]]
+    return PartitioningResult(partition, sums)
+
+
+
+
+if __name__ == "__main__":
+    import doctest
+    (failures,tests) = doctest.testmod(report=True)
+    print ("{} failures, {} tests".format(failures,tests))
diff --git a/src/numberpartitioning/solve.py b/src/numberpartitioning/solve.py
new file mode 100644
index 0000000..dd937d6
--- /dev/null
+++ b/src/numberpartitioning/solve.py
@@ -0,0 +1,115 @@
+#!python3
+
+""" 
+Utility functions for solving optimization problems using a sequence of CVXPY solvers.
+
+CVXPY supports many solvers, but some of them fail for some problems. 
+Therefore, for robustness, it may be useful to try a list of solvers, one at a time,
+   until the first one that succeeds.
+"""
+
+import cvxpy
+
+DEFAULT_SOLVERS=[cvxpy.XPRESS, cvxpy.OSQP, cvxpy.SCS]
+
+import logging
+logger = logging.getLogger(__name__)
+
+def solve(problem:cvxpy.Problem, solvers:list=DEFAULT_SOLVERS)  -> None:
+	"""
+	Try to solve the given cvxpy problem using the given solvers, in order, until one succeeds.
+    See here https://www.cvxpy.org/tutorial/advanced/index.html for a list of supported solvers.
+
+    Parameters
+    ----------
+    problem
+        a cvxpy Problem instance.
+    solvers
+        a list of one or more cvxpy solvers (optional).
+
+    Returns
+    -------
+    Nothing; the "problem" variable is automatically updated by cvxpy.
+	"""
+	is_solved=False
+	for solver in solvers[:-1]:  # Try the first n-1 solvers.
+		try:
+			problem.solve(solver=solver)
+			logger.info("Solver %s succeeds",solver)
+			is_solved = True
+			break
+		except cvxpy.SolverError as err:
+			logger.info("Solver %s fails: %s", solver, err)
+	if not is_solved:
+		problem.solve(solver=solvers[-1])   # If the first n-1 fail, try the last one.
+	if problem.status == "infeasible":
+		raise ValueError("Problem is infeasible")
+	elif problem.status == "unbounded":
+		raise ValueError("Problem is unbounded")
+
+
+def maximize(objective, constraints:list, solvers:list=DEFAULT_SOLVERS):
+	"""
+	A utility function for finding the maximum of a general objective function.
+
+    Parameters
+    ----------
+    objective
+        an expression containing cvxpy variables. Should be a concave function of the variables.
+	constraints
+	    a list of cvxpy constraints.
+    solvers
+        a list of one or more cvxpy solvers (optional).
+
+    Returns
+    -------
+	the maximum value of the objective function given the constraints.
+
+	>>> import numpy as np
+	>>> x = cvxpy.Variable()
+	>>> np.round( maximize(x, [x>=1, x<=3]), 3)
+	3.0
+	"""
+
+	problem = cvxpy.Problem(cvxpy.Maximize(objective), constraints)
+	solve(problem, solvers=solvers)
+	return objective.value.item()
+
+def minimize(objective, constraints, solvers:list=DEFAULT_SOLVERS):
+	"""
+	A utility function for finding the minimum of a general objective function.
+
+    Parameters
+    ----------
+    objective
+        an expression containing cvxpy variables. Should be a concave function of the variables.
+	constraints
+	    a list of cvxpy constraints.
+    solvers
+        a list of one or more cvxpy solvers (optional).
+
+    Returns
+    -------
+	the minimum value of the objective function given the constraints.
+
+
+	>>> import numpy as np
+	>>> x = cvxpy.Variable()
+	>>> np.round(minimize(x, [x>=1, x<=3]),3)
+	1.0
+	"""
+	problem = cvxpy.Problem(cvxpy.Minimize(objective), constraints)
+	solve(problem, solvers=solvers)
+	return objective.value.item()
+
+
+
+
+if __name__ == '__main__':
+	import sys
+	logger.addHandler(logging.StreamHandler(sys.stdout))
+	logger.setLevel(logging.INFO)
+
+	import doctest
+	(failures, tests) = doctest.testmod(report=True)
+	print("{} failures, {} tests".format(failures, tests))
diff --git a/tests/test_ilp.py b/tests/test_ilp.py
new file mode 100644
index 0000000..2ddd54b
--- /dev/null
+++ b/tests/test_ilp.py
@@ -0,0 +1,86 @@
+from numberpartitioning.ilp import *
+
+def test_max_min_partition() -> None:
+    # One part:
+    result = max_min_partition([10,20,40,0], num_parts=1, return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    # Two parts:
+    result = max_min_partition([10,20,40,0], num_parts=2, return_indices=False)
+    assert result == PartitioningResult(partition=[[10, 20, 0], [40]], sizes=[30.0, 40.0])
+
+    # Three parts:
+    result = max_min_partition([10,20,40,0], num_parts=3)
+    assert int(min(result.sizes)) == 10
+
+    # Four parts:
+    result = max_min_partition([10,20,40,0], num_parts=4)
+    assert int(min(result.sizes)) == 0
+
+    # More parts than numbers:
+    result = max_min_partition([10,20,40,0], num_parts=5)
+    assert int(min(result.sizes)) == 0
+
+    # Test multiple copies:
+    result = max_min_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    # Test different numbers of copies:
+    result = max_min_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
+
+    # Test more than one smallest-part:
+    result = max_min_partition([10,20,40,1], num_parts=3, num_smallest_parts=2)
+    assert result == PartitioningResult(partition=[[], [10, 20, 1], [40]], sizes=[0.0, 31.0, 40.0])
+
+
+
+
+def test_min_max_partition() -> None:
+    # One part:
+    result = min_max_partition([10,20,40,0], num_parts=1, return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    # Two parts:
+    result = min_max_partition([10,20,40,0], num_parts=2, return_indices=False)
+    assert result == PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])
+
+    # Three parts:
+    result = min_max_partition([10,20,40,0], num_parts=3)
+    assert int(max(result.sizes)) == 40
+
+    # More parts than numbers:
+    result = min_max_partition([10,20,40,0], num_parts=5)
+    assert int(max(result.sizes)) == 40
+
+    # Test multiple copies:
+    result = min_max_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    # Test different numbers of copies:
+    result = min_max_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
+
+    # Test more than one largest-part:
+    result = min_max_partition([10,20,40,1], num_parts=3, num_largest_parts=2)
+    PartitioningResult(partition=[[10, 1], [20], [40]], sizes=[11.0, 20.0, 40.0])
+
+
+
+
+def test_min_diff_partition() -> None:
+    # One part:
+    result = min_diff_partition([10,20,40,0], num_parts=1, return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])
+
+    # Two parts:
+    result = min_diff_partition([10,20,40,0], num_parts=2, return_indices=False)
+    assert result == PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])
+
+    # Test multiple copies:
+    result = min_diff_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
+    PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])
+
+    # Test different numbers of copies:
+    result = min_diff_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
+    assert result == PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])