diff options
Diffstat (limited to '')
| -rw-r--r-- | sumofsquares/cvxopt.py | 107 |
1 files changed, 62 insertions, 45 deletions
diff --git a/sumofsquares/cvxopt.py b/sumofsquares/cvxopt.py index 4cd803c..6a1339f 100644 --- a/sumofsquares/cvxopt.py +++ b/sumofsquares/cvxopt.py @@ -11,7 +11,7 @@ from polymatrix.expression.expression import Expression from polymatrix.polymatrix.index import MonomialIndex, VariableIndex from .abc import Problem, SolverInfo -from .constraints import NonNegative, EqualToZero +from .constraints import NonNegative, EqualToZero, PositiveSemiDefinite, ExponentialCone from .variable import OptVariable @@ -25,7 +25,6 @@ class CVXOPTInfo(SolverInfo, UserDict): return self[key] - def solve_sos_cone(prob: Problem, verbose: bool = False, *args, **kwargs) -> tuple[dict[OptVariable, float], CVXOPTInfo]: r""" Solve a conic problem in the cone of SOS polynomials :math:`\mathbf{\Sigma}_d(x)`. @@ -85,47 +84,54 @@ def solve_sos_cone(prob: Problem, verbose: bool = False, *args, **kwargs) -> tup x = poly.v_stack(prob.polynomial_variables) for c in prob.constraints: - # Convert constraints to affine expression in decision variables - constr = poly.to_affine(c.expression.quadratic_in(x).read(prob.state)) - if constr.degree > 1: - raise ValueError("CVXOpt can solve problem with constraints that are " - f"at most linear, but {str(c.expression)} has degree {constr.degree}.") - - # For each constraint we add PSD matrix - nrows, ncols = constr.shape - s_dims.append(nrows) - - # Constant term - c_coeff = constr.affine_coefficient(MonomialIndex.constant()) - c_stacked = c_coeff.T.reshape((nrows * ncols,)) - - # Linear term - l_stacked_rows = [] - for v in variable_indices.values(): - # Get the affine coefficient associated to this variable - linear = MonomialIndex((v,)) - v_coeff = constr.affine_coefficient(linear) - - # stack the columns of coeff matrix into a fat row vector (column major order). - v_stacked = v_coeff.T.reshape((nrows * ncols,)) - l_stacked_rows.append(v_stacked) - - l_stacked = np.vstack(l_stacked_rows) - - if isinstance(c, EqualToZero): - A_rows.append(l_stacked) - # Factor of -1 because it is on the right hand side - b_rows.append(-1 * c_stacked) - - elif isinstance(c, NonNegative): - if c.domain: - raise ValueError("Cannot convert non-negativity constraint to CVXOPT format. " - "Domain restricted non-negativity must be " - "preprocessed by using a Positivstellensatz!") - - # Factor of -1 because it is on the right hand side - G_rows.append(-1 * l_stacked) - h_rows.append(c_stacked) + if isinstance(c, PositiveSemiDefinite): + raise NotImplementedError + + elif isinstance(c, ExponentialCone): + raise NotImplementedError + + elif isinstance(c, EqualToZero | NonNegative): + # Convert constraints to affine expression in decision variables + constr = poly.to_affine(c.expression.quadratic_in(x).read(prob.state)) + if constr.degree > 1: + raise ValueError("CVXOpt can solve problem with constraints that are " + f"at most linear, but {str(c.expression)} has degree {constr.degree}.") + + # For each constraint we add PSD matrix + nrows, ncols = constr.shape + s_dims.append(nrows) + + # Constant term + c_coeff = constr.affine_coefficient(MonomialIndex.constant()) + c_stacked = c_coeff.T.reshape((nrows * ncols,)) + + # Linear term + l_stacked_rows = [] + for v in variable_indices.values(): + # Get the affine coefficient associated to this variable + linear = MonomialIndex((v,)) + v_coeff = constr.affine_coefficient(linear) + + # stack the columns of coeff matrix into a fat row vector (column major order). + v_stacked = v_coeff.T.reshape((nrows * ncols,)) + l_stacked_rows.append(v_stacked) + + l_stacked = np.vstack(l_stacked_rows) + + if isinstance(c, EqualToZero): + A_rows.append(l_stacked) + # Factor of -1 because it is on the right hand side + b_rows.append(-1 * c_stacked) + + elif isinstance(c, NonNegative): + if c.domain: + raise ValueError("Cannot convert non-negativity constraint to CVXOPT format. " + "Domain restricted non-negativity must be " + "preprocessed by using a Positivstellensatz!") + + # Factor of -1 because it is on the right hand side + G_rows.append(-1 * l_stacked) + h_rows.append(c_stacked) else: raise NotImplementedError(f"Cannot convert constraint of type {type(c)} " @@ -139,17 +145,28 @@ def solve_sos_cone(prob: Problem, verbose: bool = False, *args, **kwargs) -> tup b = cvxopt.matrix(np.hstack(b_rows)) if b_rows else None dims = {'l': l_dims, 'q': q_dims, 's': s_dims} + # print(q) + # print(G) + # print(h) + # print(A) + # print(b) + # print(dims) + cvxopt.solvers.options['show_progress'] = verbose if not any((G, h, A, b)): raise ValueError("Optimization problem is unconstrained!") + # Add matrices for debugging + # info = {} + info = {"q": q, "G": G, "h": h, "A": A, "b": b, "dims": dims} + if is_qp: P = cvxopt.matrix(P) - info = cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, + info |= cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, dims=dims, *args, **kwargs) else: - info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, + info |= cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, *args, **kwargs) results = {} |