diff options
| -rw-r--r-- | sumofsquares/problems.py | 44 | ||||
| -rw-r--r-- | sumofsquares/solver/cvxopt.py | 18 |
2 files changed, 40 insertions, 22 deletions
diff --git a/sumofsquares/problems.py b/sumofsquares/problems.py index 4482c1e..c313419 100644 --- a/sumofsquares/problems.py +++ b/sumofsquares/problems.py @@ -67,13 +67,15 @@ class ConicProblem(Problem): types of constraints :: - z zero cone (s = 0) - l linear cone (s >= 0) - b box cone (tl <= s <= tu) - q second order cone (|s|_2 <= t) - s positive semidefinite cone (s is PSD) - ep exponential cone ((x,y,z) in ExpCone) - ep* dual exponential cone ((u,v,w) in DualExpCone) + NAME LONG NAME MEANING + ---- --------------------------- ---------------------- + z zero cone s = 0 + l linear cone s >= 0 + b box cone tl <= s <= tu + q second order cone |s|_2 <= t + s positive semidefinite cone s is PSD + e exponential cone (x,y,z) in ExpCone + e* dual exponential cone (u,v,w) in DualExpCone p power cone p* dual power cone @@ -231,6 +233,7 @@ class SOSProblem(Problem): # Polynomial non-negativity contraints is converted to PSD # constraint of SOS quadratic form if deg.scalar().constant() > 1: + # TODO: it seems to work fine even without .symmetric(). Why? state, pm = c.expression.quadratic_in(x).symmetric().apply(state) constraints.append(PositiveSemiDefinite(pm)) @@ -239,16 +242,28 @@ class SOSProblem(Problem): state, pm = c.expression.apply(state) constraints.append(replace(c, expression=pm)) - # PSD constraint can be passed as-is elif isinstance(c, PositiveSemiDefinite): state, pm = c.expression.apply(state) + nrows, ncols = pm.shape + if nrows != ncols: + raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!") + + # PSD constraint can be passed as-is constraints.append(replace(c, expression=pm)) elif isinstance(c, ExponentialCone): - raise NotImplementedError + state, pm = c.expression.apply(state) + nrows, ncols = pm.shape + + if ncols != 3: + raise ValueError("Conic constraint must be a row vector [x, y, z] ", + "or for multiple constraints it must be an n x 3 " + f"matrix! Given expression has wrong shape {pm.shape}.") + + constraints.append(replace(c, expresssion=pm)) else: - raise NotImplementedError + raise NotImplementedError(f"Cannot process constraint of type {type(c)} (yet).") return state, InternalSOSProblem(cost, tuple(constraints), tuple(variables), polynomial_variables, @@ -291,13 +306,13 @@ class InternalSOSProblem(Problem): # See docstring of ConicProblem.constraints dims: dict[str, list[int]] = { "z": [], "l": [], "b": [], "q": [], "s": [], - "ep": [], "ep*": [], "p": [], "p*": [], + "e": [], "e*": [], "p": [], "p*": [], } # See docstring of ConicProblem.constraints constraints: dict[str, list[tuple[tuple[NDArray, ...], NDArray]]] = { "z": [], "l": [], "b": [], "q": [], "s": [], - "ep": [], "ep*": [], "p": [], "p*": [], + "e": [], "e*": [], "p": [], "p*": [], } # indices of variables in the optimization problem, sorted @@ -352,9 +367,6 @@ class InternalSOSProblem(Problem): constraints["l"].append((linear, constant)) elif isinstance(c, PositiveSemiDefinite): - if nrows != ncols: - raise ValueError(f"PSD constraint cannot contain non-square matrix of shape ({nrows, ncols})!") - dims["s"].append(nrows) constraints["s"].append((linear, constant)) @@ -362,7 +374,7 @@ class InternalSOSProblem(Problem): # Interpret row-wise dims["ep"].extend(1 for _ in range(nrows)) for i in range(nrows): - constraints["ep"].append(( + constraints["e"].append(( tuple(m[i, :] for m in linear), constant[i, :])) diff --git a/sumofsquares/solver/cvxopt.py b/sumofsquares/solver/cvxopt.py index 455ee5b..1c37fbe 100644 --- a/sumofsquares/solver/cvxopt.py +++ b/sumofsquares/solver/cvxopt.py @@ -44,10 +44,10 @@ def vectorize_matrix(m: NDArray) -> NDArray: def solve_cone(prob: ConicProblem, verbose: bool = False, *args, **kwargs) -> tuple[dict[OptVariable, NDArray | float], CVXOPTInfo]: r""" - Any `*args` and `**kwargs` other than `prob` are passed to the CVXOPT - solver. + Any `*args` and `**kwargs` other than `prob` and `vebose` are passed to the + CVXOPT solver. """ - # CVXOPT can solve problems that have the form + # CVXOPT can solve problems that have the (primal) form # # minimize .5 * x.T @ P @ x + q.T @ x # @@ -61,7 +61,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False, raise NotImplementedError("Automatic conversion of box constraints " "into linear constraints is not supported yet.") - if prob.constraints["ep"] or prob.constraints["ep*"]: + if prob.constraints["e"] or prob.constraints["e*"]: raise NotSupportedBySolver("CVXOPT cannot solve problems with the exponential cone.") if prob.constraints["p"] or prob.constraints["p*"]: @@ -121,6 +121,13 @@ def solve_cone(prob: ConicProblem, verbose: bool = False, # Solve the problem # pprint({"A": A, "b": b, "G": G, "h": h}) + # pprint(dims) + + # print(A) + # print(b) + + # print(G) + # print(h) try: if prob.is_qp: @@ -128,7 +135,6 @@ def solve_cone(prob: ConicProblem, verbose: bool = False, info = cvxopt.solvers.coneqp(P=P, q=q, G=G, h=h, A=A, b=b, dims=dims, *args, **kwargs) else: - assert q is not None, "LP has no objective! Something is very wrong" info = cvxopt.solvers.conelp(c=q, G=G, h=h, A=A, b=b, *args, **kwargs) @@ -146,7 +152,7 @@ def solve_cone(prob: ConicProblem, verbose: bool = False, results, i = {}, 0 for variable in prob.variables: num_indices = math.prod(variable.shape) - values = np.array(info['x'][i:i+num_indices]).reshape(variable.shape) + values = np.array(info["x"][i:i+num_indices]).reshape(variable.shape) if values.shape == (1, 1): values = values[0, 0] |