IntervalConstraintProgramming.jl
does not work with ModelingToolkit
v4 or higher. Use v3 instead.Download the notebook for this tutorial.
The IntervalConstraintProgramming.jl
package can be installed with
julia> using Pkg; Pkg.add("IntervalConstraintProgramming")
Once the package is installed, it can be imported. Note that you will need also the IntervalArithmetic.jl
package.
using IntervalArithmetic, IntervalConstraintProgramming
Given a domain and a set of constraint , constraint programming aims to determine all points in the domain that satisfy the constraint .
This Julia package allows you to specify a set of constraints on real-valued variables, given by (in)equalities, and rigorously calculate inner and outer approximations to the feasible set, i.e. the set that satisfies the constraints. This uses interval arithmetic provided by the IntervalArithmetic.jl
package, in particular multi-dimensional IntervalBoxes, i.e. Cartesian products of one-dimensional intervals.
The simple way to define a constraint is to use the @constraint
macro, as the following example demonstrates
S = @constraint x^2 + y^2 <= 1
Separator:
- variables: x, y
- expression: x ^ 2 + y ^ 2 ∈ [-∞, 1]
As it can be noticed, the macro itself can figure out that and are variables and you do not need to define those separately. The output of the macro is an object of type Separator
. To understand what it means, we first need to introduce a few terms.
Inner contractor: The smallest box containing all points in the domain that satisfy the constraint .
Outer contractor: The smallest box containing all points in the domain that do not satisfy the constraint .
The separator is now simply a function that returns the inner and outer contractor when applied to a domain. The domain should have always type IntervalBox
, even when it is a 1-dimensional box (i.e. an interval). Let us look an example on how to get an inner and outer contractor.
X = IntervalBox(-100..100, 2)
inner, outer = S(X)
@show inner
@show outer
inner = [-1, 1] × [-1, 1]
outer = [-100, 100] × [-100, 100]
Note Points on the boundary are counted both as inner and outer points. In other words if our constraint is e.g. , then the inner contractor will be the smallest box containing all points in the domain that satisfy and the outer contractor all points satisfying . Not that the equality is kept in both inequalities. To better understand this, consider the following example: let us look for the points in the real line satisfying , using as domain. The "pen and paper" solution would be that is the inner contractor and is the outer contractor. However, we get the following result:
S = @constraint -1<=x<=0
X = IntervalBox(-1..1, 1)
@show inner, outer = S(X)
(inner, outer) = S(X) = ([-1, 0], [-1, 1])
The inner contractor is what we expected. For the outer contractor, the key point is that the negation of is . I.e. the outer contractor looks for the smallest box inside the domain which satisfies . This is the smallest box containing and , i.e. the whole domain . Let us take another example with the same constraint but a different domain
X = IntervalBox(-0.5..0.5, 1)
@show inner, outer = S(X)
(inner, outer) = S(X) = ([-0.5, 0], [0, 0.5])
Now the outer contractor looks for the point inside satisfying , i.e .
Given a function and a set , set inversion means to find the set , which is refered as feasible set. Generally, the image set Y can be represented as a set of constraints on the expression and solving the set inversion problem means to find the points in for which these constraints are satisfied. As we show in this example, set inversion can be solved using interval constraint programming. Let's consider the function and let's find the set for which , i.e. we must find the points for which . Let's define this constraint
S = @constraint 1 <= x^2+y^2 <= 3
Separator:
- variables: x, y
- expression: x ^ 2 + y ^ 2 ∈ [1, 3]
The strategy to solve this problem is the following, we start with a large interval box, which we know will contain the feasible sets and we iteratively divide the box in smaller boxes, keeping track of which boxes are guaranteed to be inside the feasible sets, the boxes that are outside and the boxes that are on the boundary. IntervalConstraintProgramming
offers the function pave
, to do so, the signature of this function is
pave(S, X, tol)
where is our Separator, is our starting box and is a tolerance parameter. This function returns an object of type SubPaving
, which stores a list of boxes guaranteed to be into the feasible set in the attribute inner
, as well as a list of boxes on the boundary in the attribute boundary
.
X = IntervalBox(-100..100, 2) # our starting box
paving = pave(S, X, 0.125)
Paving:
- tolerance ϵ = 0.125
- inner approx. of length 164
- boundary approx. of length 164
As we can see, the function found in total 164 boxes inside the feasible set and 164 boxes at the boundary. We can now visualize how well paving approximates the feasible set
using Plots
plot(paving.inner, c="green", aspect_ratio=:equal, label="inner")
plot!(paving.boundary, c="gray", label="boundary")
The smallest the tolerance parameter, the better the approximation of the feasible set will be, as the following animation shows
tolerances = append!(collect(1:-0.1:0.1), collect(0.09:-0.01:0.01))
anim = @animate for tol in tolerances
paving = pave(S, X, tol)
plot(paving.inner, c="green", legend=false, title="tol=$tol", aspect_ratio=:equal)
plot!(paving.boundary, c="gray")
end
The @constraint
macro is very handy, but it can handle only simple expressions. If you want to handle more complicated expressions, you can use IntervalConstraintProgramming
with ModelingToolkit
. The following example shows how to import the package and define variables
using ModelingToolkit
vars = @variables x y
(x, y)
IntervalConstraintProgramming.jl
does not work with ModelingToolkit
v4 or higher. Use v3 instead.Once you have defined the variables, you can construct a separator using the Separator
constructor, which take an array of variables as first parameter and a function representing the constraint as second parameter, observe
f(x, y) = x + y < 3
S = Separator(vars, f)
Separator:
- variables: x, y
- expression: x + y == [-∞, 3]
After that, you can use the separator and the pave function as described before. As an example, let us draw the Julia logo solving a set inversion problem. We have already defined the variables, now the equations of the 3 circles are
circle1(x, y) = (x + √3)^2 + (y + 1)^2 - 9/4 < 0
circle2(x, y) = (x - √3)^2 + (y + 1)^2 - 9/4 < 0
circle3(x, y) = x^2 + (y - 2)^2 - 9/4 < 0
circle3 (generic function with 1 method)
ModelingToolkit
, then all inequalities should be strict, e.g. x+y <= 1
is not supported and should be x + y < 1
insteadNow we can define the corresponding separators.
S1 = Separator(vars, circle1)
S2 = Separator(vars, circle2)
S3 = Separator(vars, circle3)
Separator:
- variables: x, y
- expression: (x ^ 2 + (y - 2) ^ 2) - 2.25 == [-∞, 0]
Now we can use pave with diferent tolerances to plot the Julia logo
X = IntervalBox(-4..4, 2)
anim = @animate for tol in 2.0 .^ (0:-1:-6)
paving1 = pave(S1, X, tol)
paving2 = pave(S2, X, tol)
paving3 = pave(S3, X, tol)
plot(legend=false, aspect_ratio=:equal)
plot!(paving1.boundary, color=:gray, alpha=0.5)
plot!(paving2.boundary, color=:gray, alpha=0.5)
plot!(paving3.boundary, color=:gray, alpha=0.5)
plot!(paving1.inner, color=RGB(0.796, 0.235, 0.2))
plot!(paving2.inner, color=RGB(0.584, 0.345, 0.698))
plot!(paving3.inner, color=RGB(0.22, 0.596, 0.149))
end