Merge pull request #61 from TAMUparametric/revertas

Adds method for overdetermined active sets for geometric algorithm

Author: DKenefake

doc/multiobjective.rst

@@ -1,6 +1,6 @@
 from setuptools import find_packages, setup
 
-__version__ = "1.6.8"
+__version__ = "1.6.9"
 
 short_desc = (
     "Extensible Multiparametric Solver in Python"

setup.py

@@ -66,7 +66,8 @@ def manufacture_lambda(attempted, murder_list):
         return lambda x: x not in attempted and murder_list.check(x)
 
 
-def generate_reduce(candidate: tuple, murder_list=None, attempted=None, equality_set: Optional[Set[int]] = None) -> list:
+def generate_reduce(candidate: tuple, murder_list=None, attempted=None,
+                    equality_set: Optional[Set[int]] = None) -> list:
     if equality_set is None:
         equality_set = set()
 
@@ -164,6 +165,43 @@ def check(x) -> bool:
         return [[*active_set, i] for i in range(active_set[-1] + 1, num_constraints) if check([*active_set, i])]
 
 
+def find_sub_active_set(program: Union[MPLP_Program, MPQP_Program], active_set: List[int]) -> List[int]:
+    """
+    In the situation that there is an overdetermined active set we need to find the subset that is full rank and allows
+    for generating active sets that are well defined.
+
+
+    """
+
+    eq_cons = program.equality_indices
+    ineq_constraints = [i for i in active_set if i not in eq_cons]
+    kept_inequalities = []
+
+    current_rank = 0
+
+    if len(eq_cons) == 0:
+        current_rank = 0
+    else:
+        current_rank = numpy.linalg.matrix_rank(program.A[eq_cons])
+
+    for i in ineq_constraints:
+
+        trial_ineq = [*kept_inequalities, i]
+
+        A_block = numpy.block([[program.A[eq_cons]], [program.A[trial_ineq]]])
+
+        A_rank = numpy.linalg.matrix_rank(A_block)
+
+        if A_rank > current_rank:
+            kept_inequalities.append(i)
+            current_rank = A_rank
+
+        if current_rank == program.num_x():
+            return [*eq_cons, *kept_inequalities]
+
+    return [*eq_cons, *kept_inequalities]
+
+
 def get_facet_centers(A: numpy.ndarray, b: numpy.ndarray) -> List[Tuple[numpy.ndarray, numpy.ndarray, float]]:
     r"""
     This takes the polytope P, and finds all the chebychev centers and normal vectors of each facet and the radius.
@@ -192,7 +230,6 @@ def get_facet_centers(A: numpy.ndarray, b: numpy.ndarray) -> List[Tuple[numpy.nd
             if chev_ball is not None:
                 theta = chev_ball.sol[:-1]
                 radius = chev_ball.sol[-1]
-
         if theta is None:
             continue
 
@@ -254,6 +291,10 @@ def fathem_facet(center: numpy.ndarray, normal: numpy.ndarray, radius: float, pr
         # noinspection PyTypeChecker
         projected_set: List[int] = sol.active_set.tolist()
 
+        # if we have an overdetermined active set we need to find the subset that is full rank
+        if len(projected_set) > program.num_x():
+            projected_set = find_sub_active_set(program, projected_set)
+
         # test for accidental self inclusion
         if projected_set == current_active_set:
             continue

src/ppopt/mp_solvers/solver_utils.py

@@ -32,14 +32,22 @@ def process_cvxopt_solution(sol, equality_constraints, inequality_constraints, n
         lagrange[equality_constraints] = -numpy.array(sol['y']).flatten()
     lagrange[inequality_constraints] = -numpy.array(sol['z']).flatten()
 
-    non_zero_duals = numpy.where(lagrange != 0)[0]
-    active_set = numpy.array([i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
+    active = []
+
+    for i in range(num_constraints):
+        if abs(slack[i]) <= 10 ** -10 or lagrange[i] != 0:
+            active.append(i)
+
+    active = numpy.array(active)
+
+    # non_zero_duals = numpy.where(lagrange != 0)[0]
+    # active_set = numpy.array([i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
 
     if not get_duals:
         lagrange = None
 
     return SolverOutput(obj=sol['primal objective'], sol=numpy.array(sol['x']).flatten(), slack=slack,
-                        active_set=active_set,
+                        active_set=active,
                         dual=lagrange)
 
 

src/ppopt/solver_interface/cvxopt_interface.py

@@ -93,7 +93,13 @@ def solve_qp_daqp(Q: numpy.ndarray, c: Matrix, A: Matrix, b: Matrix,
 
     slack = b - A @ make_column(x_star)
 
-    non_zero_duals = numpy.where(lagrange != 0)[0]
-    active_set = numpy.array([i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
+    active = []
 
-    return SolverOutput(opt, x_star, slack.flatten(), numpy.array(active_set).astype('int64'), lagrange)
+    for i in range(num_constraints):
+        if abs(slack[i]) <= 10 ** -10 or lagrange[i] != 0:
+            active.append(i)
+
+    # non_zero_duals = numpy.where(lagrange != 0)[0]
+    # active_set = numpy.array([i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
+
+    return SolverOutput(opt, x_star, slack.flatten(), numpy.array(active).astype('int64'), lagrange)

src/ppopt/solver_interface/daqp_solver_interface.py

@@ -132,14 +132,7 @@ def solve_miqp_gurobi(Q: Matrix = None, c: Matrix = None, A: Matrix = None,
                 sol.dual = numpy.array(model.getAttr("Pi"))
 
         sol.slack = numpy.array(model.getAttr("Slack"))
-
-        if sol.dual is not None:
-            # we have a continuous problem, and we can get the duals directly
-            non_zero_duals = numpy.where(sol.dual != 0)[0]
-            sol.active_set = numpy.array(
-                [i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
-        else:
-            sol.active_set = numpy.where((A @ sol.sol.flatten() - b.flatten()) ** 2 < 10 ** -12)[0]
+        sol.active_set = numpy.where((A @ sol.sol.flatten() - b.flatten()) ** 2 < 10 ** -12)[0]
 
     return sol
 

src/ppopt/solver_interface/gurobi_solver_interface.py

@@ -68,18 +68,21 @@ def solve_qp_quadprog(Q: numpy.ndarray, c: numpy.ndarray, A: numpy.ndarray, b: n
         x_star = sol[0]
         opt = sol[1]
         duals = sol[4]
-        # active = []
+        active = []
         indexing = [*equality_constraints, *ineq]
         lagrange[indexing] = duals
         lagrange[ineq] = -lagrange[ineq]
 
         slack = b - A @ make_column(x_star)
+        for i in range(num_constraints):
+            if abs(slack[i]) <= 10 ** -10 or lagrange[i] != 0:
+                active.append(i)
 
-        non_zero_duals = numpy.where(lagrange != 0)[0]
-        active_set = numpy.array(
-            [i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
+        # non_zero_duals = numpy.where(lagrange != 0)[0]
+        # active_set = numpy.array(
+        #     [i for i in range(num_constraints) if i in non_zero_duals or i in equality_constraints])
 
-        return SolverOutput(opt, x_star, slack.flatten(), numpy.array(active_set).astype('int64'), lagrange)
+        return SolverOutput(opt, x_star, slack.flatten(), numpy.array(active).astype('int64'), lagrange)
 
     except ValueError:
         # just swallow the error as something happened Infeasibility or non-symmetry

src/ppopt/solver_interface/quad_prog_interface.py

@@ -3,7 +3,7 @@
 from src.ppopt.utils.chebyshev_ball import chebyshev_ball
 from src.ppopt.mp_solvers import mpqp_parrallel_combinatorial
 from src.ppopt.mp_solvers.solve_mpqp import mpqp_algorithm, solve_mpqp
-from tests.test_fixtures import qp_problem, simple_mpLP, portfolio_problem_analog, non_negative_least_squares
+from tests.test_fixtures import qp_problem, simple_mpLP, portfolio_problem_analog, non_negative_least_squares, over_determined_as_mplp
 
 
 def test_solve_mpqp_combinatorial(qp_problem):
@@ -84,6 +84,12 @@ def test_solve_mplp_combinatorial(simple_mpLP):
     assert solution is not None
     assert len(solution.critical_regions) == 4
 
+def tesT_solve_mplp_overdetermined_active_set(over_determined_as_mplp):
+
+    solution = solve_mpqp(over_determined_as_mplp, mpqp_algorithm.combinatorial)
+    assert solution is not None
+    assert len(solution.critical_regions) == 5
+
 
 def test_solve_missing_algorithm(qp_problem):
     try:

tests/other_tests/test_solve_mpqp.py

@@ -532,4 +532,28 @@ def pappas_multi_objective_2():
     m.add_constrs(x[i] >= 0 for i in range(1, 3))
     m.add_constrs(x[i] <= 100 for i in range(1, 3))
 
-    return m.formulate_problem()
+    return m.formulate_problem()
+
+@pytest.fixture()
+def over_determined_as_mplp():
+    """
+    The relaxation of this mpMILP was found to generate over determined active sets in some regions of the parametric
+    space. So it would cause problems with incomplete solutions.
+
+    """
+    A = numpy.array(
+        [[1, 1, 0, 0, 0], [0, 0, 1, 1, 0], [-1, 0, -1, 0, 0], [0, -1, 0, -1, -500], [-1, 0, 0, 0, 0], [0, -1, 0, 0, 0],
+         [0, 0, -1, 0, 0], [0, 0, 0, -1, 0], [0, 0, 0, 0, -1], [0, 0, 0, 0, 1]], float)
+    b = numpy.array([350, 600, 0, 0, 0, 0, 0, 0, 0, 1], float).reshape(-1, 1)
+    F = numpy.array([[0, 0], [0, 0], [-1, 0], [0, -1], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0]], float)
+    A_t = numpy.array([[1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, -1.0]], float)
+    b_t = numpy.array([[1000.0], [1000.0], [0.0], [0.0]], float)
+    H = numpy.zeros([5, 2])
+    binary_indices = [4]
+
+    Q = 2 * numpy.diag([153, 162, 162, 126, 0.1])
+    c = numpy.array([25, 25, 25, 25, 100]).reshape(-1, 1)
+
+    miqp_prog = MPMILP_Program(A, b, c, H, A_t, b_t, F, binary_indices)
+
+    return miqp_prog.generate_relaxed_problem()