diff options
| -rw-r--r-- | sumofsquares/error.py | 4 | ||||
| -rw-r--r-- | sumofsquares/problems.py | 24 |
2 files changed, 20 insertions, 8 deletions
diff --git a/sumofsquares/error.py b/sumofsquares/error.py index 34ea6d8..1f30b2a 100644 --- a/sumofsquares/error.py +++ b/sumofsquares/error.py @@ -5,3 +5,7 @@ Errors and exceptions raised by the sum of squares package. class NotSupportedBySolver(Exception): """ The chosen solver cannot solve this problem. """ + +class AlgebraicError(Exception): + """ TODO """ + diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py index 14a9cd9..3ec3f54 100644 --- a/sumofsquares/problems.py +++ b/sumofsquares/problems.py @@ -19,6 +19,7 @@ from polymatrix.variable.abc import Variable from .abc import Problem, Constraint, Solver, Result from .constraints import NonNegative +from .error import AlgebraicError from .solver.cvxopt import solve_sos_cone as cvxopt_solve_sos_cone from .utils import partition from .variable import OptVariable, from_name as opt_variable_from_name @@ -84,9 +85,6 @@ class PutinarSOSProblem(Problem): x = poly.v_stack((1,) + tuple(self.polynomial_variables)) for i, constr in enumerate(to_process): - # For each polynomial that defines the archimedean set - # we create a multiplier polynomial of the same degree - # FIXME: this is not optimal, look into how to determine the # smallest degree needed to make the positivstellensatz work. new_constr = constr.expression @@ -94,9 +92,19 @@ class PutinarSOSProblem(Problem): # To know the degree of the polynomial we need to evaluate the # expression, hence we use a monad here def make_multiplier(state): - state, p = domain_poly.degree().apply(state) - - d = p.at(0, 0).constant() + state, constr_pm = constr.expression.degree().apply(state) + state, domain_poly_pm = domain_poly.degree().apply(state) + + constr_deg = constr_pm.at(0, 0).constant() + domain_poly_deg = domain_poly_pm.at(0, 0).constant() + + # degree of multiplier + d = constr_deg - domain_poly_deg + if d < 0: + raise AlgebraicError("Cannot create Positivstellensatz multiplier because " + f"constraint polynomial has degree {constr_deg} and domain " + f"polynomial {domain_poly} has degree {domain_poly_deg}") + n = sum(math.comb(k + q - 1, k) for k in range(0, d+1)) # TODO: check that there is not a variable with this name already @@ -119,8 +127,8 @@ class PutinarSOSProblem(Problem): # Convert to polymatrix object to extract variables new_state = self.state - for constr in new_constraints: - new_state, _ = constr.expression.apply(new_state) + for c in new_constraints: + new_state, _ = c.expression.apply(new_state) def is_optvariable(v): return isinstance(v, OptVariable) |