Merge pull request #1799 from ValentinKaisermayer/vk-ineq-support

Adds support for inequalities

docs/Project.toml

@@ -3,11 +3,12 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

docs/src/systems/OptimizationSystem.md

@@ -8,9 +8,10 @@ OptimizationSystem
 
 ## Composition and Accessor Functions
 
-- `get_eqs(sys)` or `equations(sys)`: The equation to be minimized.
+- `get_op(sys)` or `objective(sys)`: The objective to be minimized.
 - `get_states(sys)` or `states(sys)`: The set of states for the optimization.
 - `get_ps(sys)` or `parameters(sys)`: The parameters for the optimization.
+- `get_constraints(sys)` or `constraints(sys)`: The constraints for the optimization.
 
 ## Transformations
 

docs/src/tutorials/optimization.md

@@ -1,26 +1,85 @@
 # Modeling Optimization Problems
 
-```julia
-using ModelingToolkit, Optimization, OptimizationOptimJL
+## Rosenbrock Function in 2D
+Let's solve the classical _Rosenbrock function_ in two dimensions.
 
-@variables x y
-@parameters a b
+First, we need to make some imports.
+```@example rosenbrock_2d
+using ModelingToolkit, Optimization, OptimizationOptimJL
+```
+Now we can define our optimization problem.
+```@example rosenbrock_2d
+@variables begin
+    x, [bounds = (-2.0, 2.0)]
+    y, [bounds = (-1.0, 3.0)]
+end
+@parameters a = 1 b = 1
 loss = (a - x)^2 + b * (y - x^2)^2
-@named sys = OptimizationSystem(loss,[x,y],[a,b])
+@named sys = OptimizationSystem(loss, [x, y], [a, b])
+```
 
+A visualization of the objective function is depicted below.
+```@eval
+using Plots
+x = -2:0.01:2
+y = -1:0.01:3
+contour(x, y, (x,y) -> (1 - x)^2 + 100 * (y - x^2)^2, fill=true, color=:viridis, ratio=:equal, xlims=(-2, 2))
+```
+
+### Explanation
+Every optimization problem consists of a set of _optimization variables_. In this case we create two variables. Additionally we assign _box constraints_ for each of them. In the next step we create two parameters for the problem with `@parameters`. While it is not needed to do this it makes it easier to `remake` the problem later with different values for the parameters. The _objective function_ is specified as well and finally everything is used to construct an `OptimizationSystem`.
+
+## Building and Solving the Optimization Problem
+Next the actual `OptimizationProblem` can be created. At this stage an initial guess `u0` for the optimization variables needs to be provided via map, using the symbols from before. Concrete values for the parameters of the system can also be provided or changed. However, if the parameters have default values assigned they are used automatically.
+```@example rosenbrock_2d
 u0 = [
-    x=>1.0
-    y=>2.0
+    x => 1.0
+    y => 2.0
 ]
 p = [
-    a => 6.0
-    b => 7.0
+    a => 1.0
+    b => 100.0
+]
+
+prob = OptimizationProblem(sys, u0, p, grad=true, hess=true)
+solve(prob, GradientDescent())
+```
+
+## Rosenbrock Function with Constraints
+```@example rosenbrock_2d_cstr
+using ModelingToolkit, Optimization, OptimizationOptimJL
+
+@variables begin
+    x, [bounds = (-2.0, 2.0)]
+    y, [bounds = (-1.0, 3.0)]
+end
+@parameters a = 1 b = 100
+loss = (a - x)^2 + b * (y - x^2)^2
+cons = [
+    x^2 + y^2 ≲ 1,
 ]
+@named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
+u0 = [
+    x => 1.0
+    y => 2.0
+]
+prob = OptimizationProblem(sys, u0, grad=true, hess=true)
+solve(prob, IPNewton())
+```
 
-prob = OptimizationProblem(sys,u0,p,grad=true,hess=true)
-solve(prob,Newton())
+A visualization of the objective function and the inequality constraint is depicted below.
+```@eval
+using Plots
+x = -2:0.01:2
+y = -1:0.01:3
+contour(x, y, (x,y) -> (1 - x)^2 + 100 * (y - x^2)^2, fill=true, color=:viridis, ratio=:equal, xlims=(-2, 2))
+contour!(x, y, (x, y) -> x^2 + y^2, levels=[1], color=:lightblue, line=4)
 ```
 
+### Explanation
+Equality and inequality constraints can be added to the `OptimizationSystem`. An equality constraint can be specified via and `Equation`, e.g., `x^2 + y^2 ~ 1`. While inequality constraints via an `Inequality`, e.g., `x^2 + y^2 ≲ 1`. The syntax is here `\lesssim` and `\gtrsim`.
+
+## Nested Systems
 Needs more text but it's super cool and auto-parallelizes and sparsifies too.
 Plus you can hierarchically nest systems to have it generate huge
 optimization problems.

src/ModelingToolkit.jl

@@ -67,7 +67,15 @@ import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
 for fun in [:toexpr]
     @eval begin
         function $fun(eq::Equation; kw...)
-            Expr(:(=), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
+            Expr(:call, :(==), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
+        end
+
+        function $fun(ineq::Inequality; kw...)
+            if ineq.relational_op == Symbolics.leq
+                Expr(:call, :(<=), $fun(ineq.lhs; kw...), $fun(ineq.rhs; kw...))
+            else
+                Expr(:call, :(>=), $fun(ineq.lhs; kw...), $fun(ineq.rhs; kw...))
+            end
         end
 
         $fun(eqs::AbstractArray; kw...) = map(eq -> $fun(eq; kw...), eqs)
@@ -86,6 +94,7 @@ abstract type AbstractTimeDependentSystem <: AbstractSystem end
 abstract type AbstractTimeIndependentSystem <: AbstractSystem end
 abstract type AbstractODESystem <: AbstractTimeDependentSystem end
 abstract type AbstractMultivariateSystem <: AbstractSystem end
+abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 
 """
 $(TYPEDSIGNATURES)
@@ -139,6 +148,7 @@ include("systems/jumps/jumpsystem.jl")
 include("systems/nonlinear/nonlinearsystem.jl")
 include("systems/nonlinear/modelingtoolkitize.jl")
 
+include("systems/optimization/constraints_system.jl")
 include("systems/optimization/optimizationsystem.jl")
 
 include("systems/pde/pdesystem.jl")
@@ -179,7 +189,7 @@ export OptimizationProblem, OptimizationProblemExpr, constraints
 export AutoModelingToolkit
 export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
-export NonlinearSystem, OptimizationSystem
+export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
@@ -192,7 +202,7 @@ export Equation, ConstrainedEquation
 export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variables, independent_variable, states, parameters, equations, controls,
-       observed, structure, full_equations
+       observed, structure, full_equations, objective
 export structural_simplify, expand_connections, linearize, linearization_function
 export DiscreteSystem, DiscreteProblem
 

src/constants.jl

@@ -32,3 +32,15 @@ macro constants(xs...)
                           xs,
                           toconstant) |> esc
 end
+
+"""
+Substitute all `@constants` in the given expression
+"""
+function subs_constants(eqs)
+    consts = collect_constants(eqs)
+    if !isempty(consts)
+        csubs = Dict(c => getdefault(c) for c in consts)
+        eqs = substitute(eqs, csubs)
+    end
+    return eqs
+end