diff options
| -rw-r--r-- | sumofsquares/__init__.py | 5 | ||||
| -rw-r--r-- | sumofsquares/cvxopt.py | 233 | ||||
| -rw-r--r-- | sumofsquares/impl/energyfunction2impl.py | 19 | ||||
| -rw-r--r-- | sumofsquares/impl/interiorconstraintimpl.py | 14 | ||||
| -rw-r--r-- | sumofsquares/impl/putinarepsilonimpl.py | 15 | ||||
| -rw-r--r-- | sumofsquares/init/initenergyfunction2.py | 110 | ||||
| -rw-r--r-- | sumofsquares/init/initinteriorconstraint.py | 28 | ||||
| -rw-r--r-- | sumofsquares/init/initputinarepsilon.py | 72 | ||||
| -rw-r--r-- | sumofsquares/mixins/energyfunction2mixin.py | 53 | ||||
| -rw-r--r-- | sumofsquares/mixins/interiorconstraintmixin.py | 21 | ||||
| -rw-r--r-- | sumofsquares/mixins/putinarepsilonmixin.py | 28 | ||||
| -rw-r--r-- | sumofsquares/mixins/sosexprmixin.py | 4 | ||||
| -rw-r--r-- | sumofsquares/mixins/sosexpropmixin.py | 9 |
13 files changed, 578 insertions, 33 deletions
diff --git a/sumofsquares/__init__.py b/sumofsquares/__init__.py index b44da40..fb40a64 100644 --- a/sumofsquares/__init__.py +++ b/sumofsquares/__init__.py @@ -2,5 +2,8 @@ from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr from sumofsquares.init.initsosexpr import init_sos_expr, init_param_expr, init_param_expr_from_reference from sumofsquares.init.initenergyfunction import init_energy_function +from sumofsquares.init.initenergyfunction2 import init_energy_function2 +from sumofsquares.init.initinteriorconstraint import init_interior_constraint +from sumofsquares.init.initputinarepsilon import init_putinar_epsilon -from sumofsquares.cvxopt import solve_cone, solve_sos_problem +from sumofsquares.cvxopt import solve_cone, solve_cone2, solve_sos_problem, solve_sos_problem2 diff --git a/sumofsquares/cvxopt.py b/sumofsquares/cvxopt.py index ed2146d..cc7d555 100644 --- a/sumofsquares/cvxopt.py +++ b/sumofsquares/cvxopt.py @@ -28,7 +28,17 @@ class CVXOptConeQPResult: -def solve_cone(state, variables, cost, s_inequality, l_inequality = tuple(), previous_result=None, print_info=False, print_matrix=False): +def solve_cone( + state, + variables, + cost, + s_inequality, + # q_inequality = tuple(), + l_inequality = tuple(), + previous_result=None, + print_info=False, + print_matrix=False, + ): state, matrix_repr = polymatrix.to_matrix_repr( (cost,) + l_inequality + s_inequality, polymatrix.v_stack(variables), @@ -135,6 +145,164 @@ def solve_cone(state, variables, cost, s_inequality, l_inequality = tuple(), pre return state, result +def solve_cone2( + state, + variables, + cost, + s_inequality = tuple(), + q_inequality = tuple(), + l_inequality = tuple(), + previous_result=None, + print_info=False, + print_matrix=False, + ): + + try: + state, matrix_repr = polymatrix.to_matrix_repr( + (cost,) + l_inequality + q_inequality + s_inequality, + polymatrix.v_stack(variables), + ).apply(state) + except: + state, matrix_repr = polymatrix.to_matrix_repr( + (cost,), + polymatrix.v_stack(variables), + ).apply(state) + + if len(l_inequality): + state, matrix_repr = polymatrix.to_matrix_repr( + l_inequality, + polymatrix.v_stack(variables), + ).apply(state) + + if len(q_inequality): + state, matrix_repr = polymatrix.to_matrix_repr( + q_inequality, + polymatrix.v_stack(variables), + ).apply(state) + + if len(s_inequality): + state, matrix_repr = polymatrix.to_matrix_repr( + s_inequality, + polymatrix.v_stack(variables), + ).apply(state) + + + + # maximum degree of cost function must be 2 + max_degree_cost = matrix_repr.data[0].get_max_degree() + assert max_degree_cost <= 2, f'{max_degree_cost=}' + + # maximum degree of constraint must be 1 + max_degree_constraints = max(data.get_max_degree() for data in matrix_repr.data[1:]) + assert max_degree_constraints == 1, f'{max_degree_constraints=}' + + cost = matrix_repr.data[0] + + # cost[1] is a 1 x n vector + q = cost[1] / 2 + + Gl = matrix_repr.data[1:1+len(l_inequality)] + dim_l = sum(G[0].shape[0] for G in Gl) + + Gq = matrix_repr.data[1+len(l_inequality):1+len(l_inequality)+len(q_inequality)] + dim_q = list(G[0].shape[0] for G in Gq) + + Gs = matrix_repr.data[1+len(l_inequality)+len(q_inequality):] + + def gen_square_matrix_dim(): + for Gs_raw in Gs: + dim = np.sqrt(Gs_raw[0].shape[0]) + + assert math.isclose(int(dim), dim), f'{dim=}' + + yield int(dim) + + dim_s = list(gen_square_matrix_dim()) + + G = np.vstack(tuple(-v[1] for v in matrix_repr.data[1:])) + h = np.vstack(tuple(v[0] for v in matrix_repr.data[1:])) + + if print_info: + print(f'number of variables: {G.shape[1]}') + print(f'{dim_l=}, {dim_q=}, {dim_s=}') + + if print_matrix: + print(f'{q=}') + print(f'{G=}') + print(f'{h=}') + + if previous_result is None: + primalstart = None + + else: + primalstart = { + 'x': cvxopt.matrix(previous_result.x), + 'y': cvxopt.matrix(previous_result.y), + 's': cvxopt.matrix(previous_result.s), + 'z': cvxopt.matrix(previous_result.z), + } + + if max_degree_cost == 1: + return_val = cvxopt.solvers.conelp( + c=cvxopt.matrix(q.reshape(-1, 1)), + G=cvxopt.matrix(G), + h=cvxopt.matrix(h), + dims={'l': dim_l, 'q': dim_q, 's': dim_s}, + primalstart=primalstart, + ) + + else: + # cost[2] is a 1 x n*n vector + nP = int(np.sqrt(cost[2].shape[1])) + P = cost[2].reshape(nP, nP).toarray() + + if print_matrix: + print(f'{P=}') + + if print_matrix: + print(f'{dim_l=}, {dim_q=}, {dim_s=}') + print(f'{P=}') + print(f'{q=}') + print(f'{G=}') + print(f'{h=}') + + return_val = cvxopt.solvers.coneqp( + P=cvxopt.matrix(P), + q=cvxopt.matrix(q.reshape(-1, 1)), + G=cvxopt.matrix(G), + h=cvxopt.matrix(h), + dims={'l': dim_l, 'q': dim_q, 's': dim_s,}, + initvals=primalstart, + ) + + x = np.array(return_val['x']).reshape(-1) + + # if x[0] == None: + # x0 = {} + # else: + # x0 = {var.name: matrix_repr.get_value(var, x) for var in variables} + + # result = CVXOptConeQPResult( + # x=x, + # y=np.array(return_val['y']).reshape(-1), + # s=np.array(return_val['s']).reshape(-1), + # z=np.array(return_val['z']).reshape(-1), + # status=return_val['status'], + # gap=return_val['gap'], + # relative_gap=return_val['relative gap'], + # primal_objective=return_val['primal objective'], + # dual_objective=return_val['dual objective'], + # primal_infeasibility=return_val['primal infeasibility'], + # dual_infeasibility=return_val['dual infeasibility'], + # primal_slack=return_val['primal slack'], + # dual_slack=return_val['dual slack'], + # iterations=return_val['iterations'], + # x0=x0 + # ) + + # return state, result + return matrix_repr.get_value(variables[0], x), return_val['status'] + def solve_sos_problem( cost: tuple[polymatrix.Expression], @@ -175,42 +343,43 @@ def solve_sos_problem( return state, result -# def solve_sos_problem( -# cost: tuple[polymatrix.Expression], -# sos_constraints: tuple[SOSConstraint], -# state: polymatrix.ExpressionState, -# subs: tuple[ParamSOSExpr] = None, -# x0: dict[ParamSOSExpr, np.ndarray] = None, -# ): -# if x0 is None: -# x0 = {} +def solve_sos_problem2( + cost: tuple[polymatrix.Expression], + sos_constraints: tuple[SOSConstraint], + state: polymatrix.ExpressionState, + subs: tuple[ParamSOSExpr] = None, + x0: dict[ParamSOSExpr, np.ndarray] = None, +): + if x0 is None: + x0 = {} -# if subs is None: -# subs = tuple() + if subs is None: + subs = tuple() -# def gen_free_param_expr(): -# for sos_constraint in sos_constraints: -# for param_expr in sos_constraint.dependence: -# if param_expr not in subs: -# yield param_expr + def gen_free_param_expr(): + for sos_constraint in sos_constraints: + for param_expr in sos_constraint.dependence: + if param_expr not in subs: + yield param_expr -# free_param_expr = tuple(set(gen_free_param_expr())) + free_param_expr = tuple(set(gen_free_param_expr())) + + # print(f'{tuple(val.param.name for val in free_param_expr)=}') -# sub_vals = tuple((sub.param, x0[sub.param]) for sub in subs) + sub_vals = tuple((sub.param, x0[sub.param.name]) for sub in subs if sub is not None and sub.param.name in x0) -# def gen_inequality(): -# for sos_constraint in sos_constraints: -# if any(param_expr in free_param_expr for param_expr in sos_constraint.dependence): -# # print(tuple(d.name for d in sos_constraint.dependence)) -# yield sos_constraint.constraint.eval(sub_vals) + def gen_inequality(): + for sos_constraint in sos_constraints: + if any(param_expr in free_param_expr for param_expr in sos_constraint.dependence): + yield sos_constraint.constraint.eval(sub_vals) -# inequality = tuple(gen_inequality()) + inequality = tuple(gen_inequality()) -# state, result = solve_cone( -# state, -# variables=tuple(param_expr.param for param_expr in free_param_expr), -# cost=sum(cost), -# s_inequality=inequality, -# ) + state, result = solve_cone( + state, + variables=tuple(param_expr.param for param_expr in free_param_expr), + cost=sum(cost), + s_inequality=inequality, + ) -# return state, result + return state, result diff --git a/sumofsquares/impl/energyfunction2impl.py b/sumofsquares/impl/energyfunction2impl.py new file mode 100644 index 0000000..7e6b091 --- /dev/null +++ b/sumofsquares/impl/energyfunction2impl.py @@ -0,0 +1,19 @@ +import dataclassabc +import polymatrix + +from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr +from sumofsquares.mixins.energyfunction2mixin import EnergyFunction2Mixin +from sumofsquares.sosconstraint.abc import SOSConstraint + + +@dataclassabc.dataclassabc(frozen=True) +class EnergyFunction2Impl(EnergyFunction2Mixin): + name: str + h: SOSExpr + h_dot: SOSExpr + epsilon: ParamSOSExpr + gamma_b: ParamSOSExpr + gamma_roi: list[ParamSOSExpr] + # gamma_w_inner: ParamSOSExpr | None + # gamma_w_outer: ParamSOSExpr | None + sos_constraints: tuple[SOSConstraint] diff --git a/sumofsquares/impl/interiorconstraintimpl.py b/sumofsquares/impl/interiorconstraintimpl.py new file mode 100644 index 0000000..ac76c19 --- /dev/null +++ b/sumofsquares/impl/interiorconstraintimpl.py @@ -0,0 +1,14 @@ +import dataclassabc +import polymatrix + +from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr +from sumofsquares.mixins.energyfunction2mixin import EnergyFunction2Mixin +from sumofsquares.mixins.interiorconstraintmixin import InteriorConstraintMixin +from sumofsquares.sosconstraint.abc import SOSConstraint + + +@dataclassabc.dataclassabc(frozen=True) +class InteriorConstraintImpl(InteriorConstraintMixin): + name: str + gamma: ParamSOSExpr | None + sos_constraints: tuple[SOSConstraint] diff --git a/sumofsquares/impl/putinarepsilonimpl.py b/sumofsquares/impl/putinarepsilonimpl.py new file mode 100644 index 0000000..dfba3f0 --- /dev/null +++ b/sumofsquares/impl/putinarepsilonimpl.py @@ -0,0 +1,15 @@ +import dataclassabc +import polymatrix + +from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr +from sumofsquares.mixins.putinarepsilonmixin import PutinarEpsilonMixin +from sumofsquares.sosconstraint.abc import SOSConstraint + + +@dataclassabc.dataclassabc(frozen=True) +class PutinarEpsilonImpl(PutinarEpsilonMixin): + name: str + epsilon: ParamSOSExpr + gamma: dict[str, ParamSOSExpr] + sos_constraints: tuple[SOSConstraint] + condition: SOSExpr diff --git a/sumofsquares/init/initenergyfunction2.py b/sumofsquares/init/initenergyfunction2.py new file mode 100644 index 0000000..785ef87 --- /dev/null +++ b/sumofsquares/init/initenergyfunction2.py @@ -0,0 +1,110 @@ +import polymatrix + +from sumofsquares.abc.sosexpr import SOSExpr +from sumofsquares.impl.energyfunction2impl import EnergyFunction2Impl +from sumofsquares.init.initsosexpr import init_param_expr, init_param_expr_from_reference + + +def init_energy_function2( + name: str, + # h: SOSExpr, + # x_dot: SOSExpr, + cond: SOSExpr, + h: SOSExpr, + positive_b: bool, + epsilon_min: SOSExpr | None = None, + roi: list[SOSExpr | polymatrix.Expression] = None, + decrease_rate: SOSExpr | polymatrix.Expression | None = None, + # w_outer: SOSExpr | polymatrix.Expression | None = None, + # w_inner: SOSExpr | polymatrix.Expression | None = None, +): + sos_constraints = tuple() + + cond = cond.cache() + + # nh = h.diff().T + # h_dot = (nh.T @ x_dot).cache() + + gamma_b = init_param_expr_from_reference( + name=f'gamma_b_{name}', + reference=cond, + multiplicand=h, + ) + if positive_b: + sos_constraints += gamma_b.sos_constraints + + energy_expr = -cond + energy_expr -= gamma_b * h + + if epsilon_min is None: + epsilon = None + + else: + epsilon = init_param_expr( + name=f'epsilon_{name}', + variables=h.variables, + ) + + sos_constraints += (epsilon - epsilon_min).sos_constraints + energy_expr += epsilon + + if roi is None: + gamma_roi = None, + + else: + def gen_gamma_roi(): + for idx, roi_i in enumerate(roi): + yield init_param_expr_from_reference( + name=f'gamma_roi{idx}_{name}', + reference=cond, + multiplicand=roi_i, + ) + gamma_roi = tuple(gen_gamma_roi()) + + for gamma_roi_i, roi_i in zip(gamma_roi, roi): + sos_constraints += gamma_roi_i.sos_constraints + + energy_expr += gamma_roi_i * roi_i + + if decrease_rate is not None: + energy_expr -= decrease_rate + + sos_constraints += energy_expr.sos_constraints + + # if w_outer is None: + # gamma_w_outer = None + + # else: + # gamma_w_outer = init_param_expr_from_reference( + # name=f'gamma_w_outer_{name}', + # reference=h, + # multiplicand=w_outer, + # ) + + # sos_constraints += gamma_w_outer.sos_constraints + # sos_constraints += (h - gamma_w_outer * w_outer).sos_constraints + + # if w_inner is None: + # gamma_w_inner = None + + # else: + # gamma_w_inner = init_param_expr_from_reference( + # name=f'gamma_w_inner_{name}', + # reference=h, + # multiplicand=w_inner, + # ) + + # sos_constraints += gamma_w_inner.sos_constraints + # sos_constraints += (gamma_w_inner * w_inner - h).sos_constraints + + return EnergyFunction2Impl( + name=name, + h=h, + h_dot=cond, + epsilon=epsilon, + gamma_roi=gamma_roi, + gamma_b=gamma_b, + # gamma_w_outer=gamma_w_outer, + # gamma_w_inner=gamma_w_inner, + sos_constraints=sos_constraints, + ) diff --git a/sumofsquares/init/initinteriorconstraint.py b/sumofsquares/init/initinteriorconstraint.py new file mode 100644 index 0000000..3daa44e --- /dev/null +++ b/sumofsquares/init/initinteriorconstraint.py @@ -0,0 +1,28 @@ +import polymatrix + +from sumofsquares.abc.sosexpr import SOSExpr +from sumofsquares.impl.interiorconstraintimpl import InteriorConstraintImpl +from sumofsquares.init.initsosexpr import init_param_expr_from_reference + + +def init_interior_constraint( + name: str, + h: SOSExpr, + w: SOSExpr | polymatrix.Expression, +): + sos_constraints = tuple() + + gamma = init_param_expr_from_reference( + name=f'gamma_{name}', + reference=h, + multiplicand=w, + ) + + sos_constraints += gamma.sos_constraints + sos_constraints += (h - gamma * w).sos_constraints + + return InteriorConstraintImpl( + name=name, + gamma=gamma, + sos_constraints=sos_constraints, + ) diff --git a/sumofsquares/init/initputinarepsilon.py b/sumofsquares/init/initputinarepsilon.py new file mode 100644 index 0000000..d3e348d --- /dev/null +++ b/sumofsquares/init/initputinarepsilon.py @@ -0,0 +1,72 @@ +import polymatrix + +from sumofsquares.abc.sosexpr import SOSExpr +from sumofsquares.impl.energyfunction2impl import EnergyFunction2Impl +from sumofsquares.impl.putinarepsilonimpl import PutinarEpsilonImpl +from sumofsquares.init.initsosexpr import init_param_expr, init_param_expr_from_reference + + +def init_putinar_epsilon( + name: str, + f0: SOSExpr, + fi: dict[str, SOSExpr | polymatrix.Expression] = None, + gi: dict[str, SOSExpr | polymatrix.Expression] = None, + epsilon_min: SOSExpr | None = None, + decrease_rate: SOSExpr | polymatrix.Expression | None = None, +): + sos_constraints = tuple() + + condition = f0 + + def gen_gamma(fi): + for key, val in fi.items(): + yield key, init_param_expr_from_reference( + name=f'gamma_{key}_{name}', + reference=f0, + # reference=condition, + multiplicand=val, + ) + + if fi is None: + gamma_fi = {} + + else: + gamma_fi = dict(gen_gamma(fi)) + + for key, gamma_i in gamma_fi.items(): + sos_constraints += gamma_i.sos_constraints + condition = condition - gamma_i * fi[key] + + if gi is None: + gamma_gi = {} + + else: + gamma_gi = dict(gen_gamma(gi)) + + for key, gamma_i in gamma_gi.items(): + condition = condition - gamma_i * gi[key] + + if epsilon_min is None: + epsilon = None + + else: + epsilon = init_param_expr( + name=f'epsilon_{name}', + variables=f0.variables, + ) + + sos_constraints += (epsilon - epsilon_min).sos_constraints + condition += epsilon + + if decrease_rate is not None: + condition -= decrease_rate + + sos_constraints += condition.sos_constraints + + return PutinarEpsilonImpl( + name=name, + epsilon=epsilon, + gamma=gamma_fi | gamma_gi, + sos_constraints=sos_constraints, + condition=condition, + ) diff --git a/sumofsquares/mixins/energyfunction2mixin.py b/sumofsquares/mixins/energyfunction2mixin.py new file mode 100644 index 0000000..f96abb6 --- /dev/null +++ b/sumofsquares/mixins/energyfunction2mixin.py @@ -0,0 +1,53 @@ +import abc +import polymatrix + +from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr +from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin +from sumofsquares.sosconstraint.abc import SOSConstraint + + +class EnergyFunction2Mixin(GetSOSConstraintMixin, abc.ABC): + @property + @abc.abstractmethod + def name(self) -> str: + ... + + @property + @abc.abstractmethod + def h(self) -> SOSExpr: + ... + + @property + @abc.abstractmethod + def h_dot(self) -> SOSExpr: + ... + + # @property + # @abc.abstractmethod + # def variables(self) -> polymatrix.Expression: + # ... + + @property + @abc.abstractmethod + def epsilon(self) -> ParamSOSExpr: + ... + + @property + @abc.abstractmethod + def gamma_b(self) -> ParamSOSExpr: + ... + + @property + @abc.abstractmethod + def gamma_roi(self) -> list[ParamSOSExpr]: + ... + + # @property + # @abc.abstractmethod + # def gamma_w_outer(self) -> ParamSOSExpr | None: + # ... + + # @property + # @abc.abstractmethod + # def gamma_w_inner(self) -> ParamSOSExpr | None: + # ... diff --git a/sumofsquares/mixins/interiorconstraintmixin.py b/sumofsquares/mixins/interiorconstraintmixin.py new file mode 100644 index 0000000..781a3d4 --- /dev/null +++ b/sumofsquares/mixins/interiorconstraintmixin.py @@ -0,0 +1,21 @@ +import abc + +from sumofsquares.abc.sosexpr import ParamSOSExpr +from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin + + +class InteriorConstraintMixin(GetSOSConstraintMixin, abc.ABC): + @property + @abc.abstractmethod + def name(self) -> str: + ... + + # @property + # @abc.abstractmethod + # def h(self) -> SOSExpr: + # ... + + @property + @abc.abstractmethod + def gamma(self) -> ParamSOSExpr | None: + ... diff --git a/sumofsquares/mixins/putinarepsilonmixin.py b/sumofsquares/mixins/putinarepsilonmixin.py new file mode 100644 index 0000000..6c1cd53 --- /dev/null +++ b/sumofsquares/mixins/putinarepsilonmixin.py @@ -0,0 +1,28 @@ +import abc +import polymatrix + +from sumofsquares.abc.sosexpr import ParamSOSExpr, SOSExpr +from sumofsquares.mixins.getsosconstraintmixin import GetSOSConstraintMixin +from sumofsquares.sosconstraint.abc import SOSConstraint + + +class PutinarEpsilonMixin(GetSOSConstraintMixin, abc.ABC): + @property + @abc.abstractmethod + def name(self) -> str: + ... + + @property + @abc.abstractmethod + def epsilon(self) -> ParamSOSExpr: + ... + + @property + @abc.abstractmethod + def gamma(self) -> dict[str, ParamSOSExpr]: + ... + + @property + @abc.abstractmethod + def condition(self) -> SOSExpr: + ... diff --git a/sumofsquares/mixins/sosexprmixin.py b/sumofsquares/mixins/sosexprmixin.py index d0c0d17..f3d4dc9 100644 --- a/sumofsquares/mixins/sosexprmixin.py +++ b/sumofsquares/mixins/sosexprmixin.py @@ -30,6 +30,10 @@ class SOSExprMixin(GetSOSConstraintMixin, abc.ABC): return self.underlying.sos_matrix @property + def sos_matrix_vec(self) -> polymatrix.Expression: + return self.underlying.sos_matrix_as_vector + + @property def sos_constraints(self) -> tuple[SOSConstraint]: return (init_sos_constraint( dependence=self.dependence, diff --git a/sumofsquares/mixins/sosexpropmixin.py b/sumofsquares/mixins/sosexpropmixin.py index 2020043..f7cf0c2 100644 --- a/sumofsquares/mixins/sosexpropmixin.py +++ b/sumofsquares/mixins/sosexpropmixin.py @@ -47,6 +47,10 @@ class SOSExprOPMixin(SOSExprMixin): ) else: + # var_set_left = set(left.variables) + # var_set_right = set(right.variables) + # assert var_set_left.issubset(var_set_right) or var_set_right.issubset(var_set_left), f'{left.variables=}, {right.variables=}' + assert left.variables == right.variables, f'{left.variables=}, {right.variables=}' underlying=init_sos_expr_base( @@ -142,3 +146,8 @@ class SOSExprOPMixin(SOSExprMixin): ), ) + def v_stack(self, other): + def op(left, right): + return polymatrix.v_stack((left, right)) + + return self._binary(op, self, other) |