SetProg

SetProg is a JuMP extension for Set Programming : optimization with set variables and inclusion/containment constraints. This package allows the formulation of a mathematical programming involving both classical variables and constraints supported by JuMP and set variables and constraints.

Variables on fixed domains

The Domain wrapper restricts a set variable to a fixed polyhedral region. For an orthant whose dimension is inferred from the model, use MOI.Nonnegatives:

using SetProg
model = Model()
@variable(model, S, Domain(Polytope(symmetric = true), MOI.Nonnegatives))
orthant = polyhedron(hrep([-1.0 0.0; 0.0 -1.0], zeros(2)))
unit = SetProg.Sets.Intersection(SetProg.Sets.PolarPoint(ones(2)), orthant)
@constraint(model, S ⊆ unit)

Here the unit simplex supplies dimension two and bounds the linear gauge's coefficients below by one. MOI.Nonnegatives(n) specifies the dimension explicitly. A fixed Polyhedra representation or polyhedron is also accepted as the second argument, with its dimension checked against the wrapped family. The wrapper also accepts other families, such as Domain(Ellipsoid(symmetric = true), MOI.Nonnegatives), and constructs a Sets.Intersection. Available constraints depend on the underlying set type. The constraints select the modeling space. S ⊆ unit uses the restricted primal gauge. A normalization written as polar(unit) ⊆ S selects dual modeling, returning the polar of the restricted gauge. For a positive simplex unit, this is the set of points with coordinate bounds xᵢ ≤ cᵢ, whose intersection with the nonnegative orthant is a positive box.

Sublevel sets without a positivity constraint

Use psd = false to omit the PSD constraint on an ellipsoid's coefficient matrix, or the SOS constraint on a polynomial's Gram matrix. For example, Ellipsoid(symmetric = false, psd = false, dimension = 2) represents x'Qx + 2b'x + β ≤ 0, with free symmetric Q, vector b, and scalar β. PolySet(degree = 0:4, psd = false, variables = x) represents p(x) ≤ 1 with monomials of degrees zero through four and a free Gram matrix. An integer, such as degree = 4, selects homogeneous polynomials. Degree ranges and collections work with both values of psd; their maximum must be even. symmetric = true excludes odd-degree terms. These affine families use primal modeling; homogeneous families keep the existing space selection. Defaults retain their positivity constraints.

Separation certificates

Quadratic and polynomial sublevel sets support S-procedure certificates:

using SetProg, LinearAlgebra
model = Model()
@variable(model, S, Ellipsoid(psd = false, dimension = 2))
initial = SetProg.Sets.ShiftedEllipsoid(
    LinearAlgebra.Symmetric([1.0 0.0; 0.0 1.0]), zeros(2), -16.0,
)
unsafe = SetProg.Sets.ShiftedEllipsoid(
    LinearAlgebra.Symmetric([-1.0 0.0; 0.0 -1.0]), zeros(2), 36.0,
)
@variable(model, ε >= 0)
@constraint(model, initial ⊆ S, margin = ε)
@constraint(model, [S, unsafe] in Disjoint(), margin = ε)

With zero-level polynomials a and b, inclusion certifies λa - b ≥ 0 and disjointness certifies λa + b ≥ 0 when the first set is fixed. When the second set is fixed, the multiplier scales b instead. A positive margin subtracts margin * (1 + ‖x‖²)^d from this certificate, where 2d is its degree. Quadratic certificates are PSD matrix constraints; higher degrees use SOS constraints and require using SumOfSquares.

By default a fixed set is scaled by a nonnegative scalar. Pass S_procedure_scaling = SProcedureMultiplier(degree = 2) for an SOS multiplier, or a nonnegative number for a fixed multiplier. Two variable sets require a fixed multiplier, defaulting to one. The Disjoint constraint defaults to a margin of one, certifying strict disjointness. A zero-margin certificate alone can allow the sets to share boundary points.

Piecewise-linear inclusions

For convex piecewise-linear gauges with a bounded source partition, inclusion is checked once at each vertex of that partition. These vertices generate the source cones: linearity on each source cone and convexity of the target extend the inequalities to the whole space, even when source cones cross target cones. No cone-intersection multipliers are introduced.

SetProg's polytope families impose continuity and convexity. Fixed numeric pieces are checked for convexity before using this certificate. Sets without that guarantee, or with unbounded source partitions, retain the general cone-domain certificate.

Polyhedral separation

A polytope family with a fixed partition also supports ellipsoid containment and separation from an ellipsoid's exterior:

@variable(model, S, Polytope(symmetric = true, piecewise = partition))
@constraint(model, initial ⊆ S, margin = ε)
@constraint(model, [S, unsafe] in Disjoint(), margin = ε)

Here initial is a bounded ellipsoid and unsafe is a ShiftedEllipsoid whose quadratic part is negative definite: it describes the exterior of a bounded ellipsoid. The partition must be bounded and contain the origin in its interior. Both modeling spaces are supported, and other constraints can select between them: in primal gauge coordinates, initial containment uses second-order cones and unsafe separation is linear on the fixed cone rays; in dual vertex coordinates, initial containment is linear and each vertex is constrained by an SOC. No SOS dependency is needed.

Margins are representation-dependent. In primal coordinates, the gauge is at most 1 - ε on the initial set and exceeds the safe ellipsoid's gauge by ε * norm(direction) on each cone ray. In dual coordinates, the margin pads the initial support and reduces the safe ellipsoid's Cholesky radius. A positive margin certifies strict separation in either representation. certificate_matrix assembles all conic parts for numerical checks; dual returns one dual vector per part when a constraint generates several cones.

Quadratic separation on conic partitions

A piecewise ellipsoid supports initial-set inclusion and disjointness from a quadratic sublevel set, including the exterior of an ellipsoid:

using SetProg, SumOfSquares
@variable(model, S, Ellipsoid(symmetric = true, piecewise = partition))
@constraint(model, initial ⊆ S, margin = ε)
@constraint(model, [S, unsafe] in Disjoint(), margin = ε)

The initial and unsafe sets are fixed Ellipsoid or ShiftedEllipsoid objects. Each piece uses an S-procedure certificate restricted to its cone. After homogenization, nonnegative products of the cone inequalities and a PSD remainder certify the quadratic inequality. These are sufficient SDP certificates; they do not require each quadratic piece to satisfy the inequality outside its cone. margin multiplies 1 + dot(x, x).

The default uses an independent nonnegative scalar multiplier for each piece. S_procedure_scaling can also be a vector with one multiplier per piece, or a single multiplier shared by all pieces. certificate_matrix returns the block-diagonal matrix of the fitted PSD remainders, and dual returns one dual per piece.

The same domain certificates apply directly to Sets.Intersection(quadratic_set, cone) as the subset in an inclusion or as either operand of Disjoint(). The domain must be a homogeneous polyhedral cone. This also supports shifted quadratic inequalities on cones.

Lₚ sums and polar Lₚ sums

LpSum(C, S, p) is the Firey sum, defined by its support function: h(LpSum(C, S, p), x)^p = h(C, x)^p + h(S, x)^p. PolarLpSum(C, S, p) is its polar counterpart: polar(LpSum(polar(C), polar(S), p)). Its gauge satisfies g(PolarLpSum(C, S, p), x)^p = g(C, x)^p + g(S, x)^p. The exponent must be finite and at least one. The operands must represent convex sets containing the origin with nonnegative gauges.

The exponent belongs to the operation, independently of the set family. For example, exponent two adds squared polyhedral gauges or the quadratic gauges of ellipsoids; exponent four adds the defining polynomials of homogeneous quartic PolySets. Exponent one adds ordinary gauges, while LpSum(C, S, 1) is the Minkowski sum.

using SetProg, LinearAlgebra
C = SetProg.Sets.Ellipsoid(Symmetric([1.0 0.0; 0.0 1.0]))
S = SetProg.Sets.Ellipsoid(Symmetric([2.0 0.0; 0.0 2.0]))
combined = PolarLpSum(C, S, 2)
SetProg.sublevel_eval(combined, [1.0, 2.0]) # 15.0
support_function([1.0, 2.0], LpSum(C, S, 2)) # √7.5

sublevel_eval for PolarLpSum evaluates the sum of the powered gauges; the represented set is where this sum is at most one. support_function for LpSum evaluates its ordinary, degree-one support function. Evaluating polynomial gauges requires loading SumOfSquares and homogeneous operands.

These expressions support dimension inference, polynomial coordinates, set-variable loading, and JuMP.value.

For exponent two, an ellipsoid can be constrained inside a PolarLpSum of fixed ellipsoids and polars of matrix ellipsoids, including their linear preimages. SetProg writes a Schur-complement PSD constraint, keeping inverse matrices as diagonal blocks. This supports affine maps and affine matrix ellipsoids without multiplying optimization variables. For example:

model = Model()
@variable(model, X[1:2, 1:2], Symmetric)
@variable(model, D[1:2, 1:2], Symmetric)
@variable(model, A[1:2, 1:2])
source = SetProg.Sets.Ellipsoid(X)
destination = SetProg.Sets.LinearPreImage(polar(SetProg.Sets.Ellipsoid(D)), A)
@constraint(model, source ⊆ PolarLpSum(destination, C, 2))

The optional psd_margin keyword subtracts a uniform diagonal margin from this certificate matrix; it defaults to zero. Homogeneous polynomial gauges use the SOS certificates described below.

See Section 5.5 of Gardner, Hug and Weil, Operations between sets in geometry. Their polar-addition convention uses parameter -p; these constructors take positive p for both operations.

Polynomial powered-gauge certificates

For homogeneous polynomial gauges, SetProg certifies a PolarLpSum inclusion by requiring the difference of powered gauges to be SOS. For example, exponent four adds a quartic value polynomial and the square of a quadratic cost:

using SetProg, SumOfSquares, LinearAlgebra
model = Model()
@variable(model, V, PolySet(symmetric = true, convex = true,
                          psd = false, degree = 4, dimension = 2))
F = 0.5 * Matrix{Float64}(I, 2, 2)
cost = SetProg.Sets.Ellipsoid(Symmetric(Matrix{Float64}(I, 2, 2)))
next_value = SetProg.Sets.LinearPreImage(V, F)
@constraint(model, V ⊆ PolarLpSum(next_value, cost, 4), margin = 1e-6)
@objective(model, Max, L1_heuristic(volume(V), ones(2)))

This imposes V(x) ≥ V(Fx) + (x'x)^2 + 1e-6 * (x'x)^2 and minimizes the integral of V over the unit box. The convex = true family supplies an SOS-convexity certificate. Homogeneity and convexity make V(x)^(1/4) a gauge, including when a separate PSD constraint on its defining Gram matrix is omitted with psd = false.

The polynomial certificate requires integer powers of the defining polynomials. The exponent is shared by nested sums. Nonpolynomial cases, such as exponent two applied to a quartic defining polynomial, require other certificates. The inclusion is an SDP when the residual coefficients are affine in the optimization variables; maps containing feedback variables can make those coefficients nonlinear.

Piecewise quadratic powered-gauge certificates

For Ellipsoid(symmetric = true, piecewise = partition), inclusions into PolarLpSum(..., 2) use quadratic positivity certificates on each common cone. The source cones are intersected with the preimages of the cones of every piecewise operand. On each intersection the squared gauges are quadratic forms, so their difference is constrained to be nonnegative using the same cone multipliers as piecewise quadratic inclusion.

Operands may be ellipsoids, piecewise ellipsoids, fixed linear preimages, or nested sums with exponent two. The resulting constraints are an SDP when the quadratic matrices are affine and the maps are fixed. The cone positivity certificates can be conservative. margin adds margin * ‖x‖² to the target squared gauge.

Polyhedral powered-gauge certificates

For a Polytope family with a fixed piecewise partition, inclusions into PolarLpSum use the extreme rays of each source cone. On that cone the source gauge is linear and the target gauge is sublinear, so checking the rays suffices. This preserves the supplied partition without intersecting it with the target's cones.

Exponent one produces linear constraints; exponent two produces second-order cone constraints. Operands may be piecewise-linear gauges or fixed ellipsoids, with fixed linear preimages and nested sums of the same exponent. For proper gauges these ray tests are exact. For an unbounded piecewise-linear target, the squared certificate uses the absolute value of each selected linear piece and is sufficient, potentially conservative. The optional margin adds margin * ‖x‖^p to the powered target gauge.

model = Model()
@variable(model, V, Polytope(symmetric = true, piecewise = partition))
next_value = SetProg.Sets.LinearPreImage(V, F)
@constraint(model, V ⊆ PolarLpSum(next_value, cost, 2))
@objective(model, Max, L1_heuristic(volume(V)))

Here partition is a bounded polytope containing zero in its interior, F is fixed, and cost is a fixed gauge set. The default volume heuristic minimizes integrals of the linear gauge pieces over bounded sections of their cones. Neither exponent requires SumOfSquares.