From 58f37a9d5d27d875846c4d9fdb9f9356cd99ec3e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Tue, 5 Nov 2024 13:00:27 +0100
Subject: [PATCH] Use Delta

---
 docs/src/examples/solvers/Unicycle robot.jl   | 24 ++------
 src/MOI_wrapper.jl                            | 59 +++++++++++++++----
 .../abstraction/uniform_grid_abstraction.jl   | 41 ++++++-------
 src/system/controlsystem.jl                   |  2 +-
 4 files changed, 73 insertions(+), 53 deletions(-)

diff --git a/docs/src/examples/solvers/Unicycle robot.jl b/docs/src/examples/solvers/Unicycle robot.jl
index c1aabb706..73d0efadb 100644
--- a/docs/src/examples/solvers/Unicycle robot.jl	
+++ b/docs/src/examples/solvers/Unicycle robot.jl	
@@ -48,9 +48,10 @@ x_upp = -x_low
 @variable(model, -1 <= u[1:2] <= 1)
 
 # Define the dynamics
-@constraint(model, ∂(x[1]) == x[1] + u[1] * cos(x[3]))
-@constraint(model, ∂(x[2]) == x[2] + u[1] * sin(x[3]))
-@constraint(model, ∂(x[3]) == x[3] + u[2])
+@constraint(model, Δ(x[1]) == x[1] + u[1] * cos(x[3]))
+@constraint(model, Δ(x[2]) == x[2] + u[1] * sin(x[3]))
+#@constraint(model, Δ(x[3]) == rem(x[3] + u[2], 2 * π))
+@constraint(model, Δ(x[3]) == x[3] + u[2])
 
 # Define the initial and target sets
 x_initial = [1.0, -1.7, 0.0]
@@ -129,23 +130,6 @@ end
 
 # ### Definition of the abstraction
 
-# First we need to set the mode of the abstraction to `DIONYSOS.Optim.Abstraction.UniformGridAbstraction.DSICRETE_TIME`:
-set_attribute(
-    model,
-    "approx_mode",
-    Dionysos.Optim.Abstraction.UniformGridAbstraction.DISCRETE_TIME,
-)
-
-# We define the system map $f$:
-function sys_map(x, u, _)
-    return StaticArrays.SVector{3}(
-        x[1] + u[1] * cos(x[3]),
-        x[2] + u[1] * sin(x[3]),
-        (x[3] + u[2]) % (2 * π),
-    )
-end
-set_attribute(model, "system_map", sys_map)
-
 # We define the growth bound function of $f$:
 function growth_bound(r, u, _)
     β = u[1] * r[3]
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index 31ec06bbd..f1881cf49 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -7,10 +7,13 @@ import Symbolics
 
 @enum(VariableType, INPUT, STATE, MODE)
 
+@enum(TimeType, UNKNOWN, CONTINUOUS, DISCRETE)
+
 mutable struct Optimizer <: MOI.AbstractOptimizer
     inner::Any
     nlp_model::Any
     evaluator::Union{Nothing, MOI.Nonlinear.Evaluator}
+    time_type::TimeType
     variable_values::Vector{Float64}
     derivative_values::Vector{Float64}
     variable_types::Vector{VariableType}
@@ -29,6 +32,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
             MOI.instantiate(Dionysos.Optim.Abstraction.UniformGridAbstraction.Optimizer),
             nothing,
             nothing,
+            UNKNOWN,
             Float64[],
             Float64[],
             VariableType[],
@@ -49,6 +53,7 @@ end
 MOI.is_empty(model::Optimizer) = isempty(model.variable_types)
 
 function MOI.empty!(model::Optimizer)
+    model.time_type = UNKNOWN
     empty!(model.variable_values)
     empty!(model.derivative_values)
     empty!(model.variable_types)
@@ -168,7 +173,21 @@ function MOI.add_constraint(
         if lhs isa MOI.ScalarNonlinearFunction && length(lhs.args) == 1
             var = lhs.args[]
             if var isa MOI.VariableIndex
-                if lhs.head == :diff
+                if lhs.head == :∂
+                    if model.time_type == DISCRETE
+                        error("Cannot add constraint with `∂` since you already added a constraint with `Δ`.")
+                    end
+                    model.time_type = CONTINUOUS
+                    model.dynamic[var.value] = rhs
+                    model.nonlinear_index += 1
+                    return MOI.ConstraintIndex{typeof(func), typeof(set)}(
+                        model.nonlinear_index,
+                    )
+                elseif lhs.head == :Δ
+                    if model.time_type == CONTINUOUS
+                        error("Cannot add constraint with `Δ` since you already added a constraint with `∂`.")
+                    end
+                    model.time_type = DISCRETE
                     model.dynamic[var.value] = rhs
                     model.nonlinear_index += 1
                     return MOI.ConstraintIndex{typeof(func), typeof(set)}(
@@ -303,13 +322,23 @@ function system(
     )
     _U_ =
         Dionysos.Utils.HyperRectangle(_svec(model.lower, u_idx), _svec(model.upper, u_idx))
-    return MathematicalSystems.ConstrainedBlackBoxControlContinuousSystem(
-        dynamic(model, x_idx, u_idx),
-        Dionysos.Utils.get_dims(_X_),
-        Dionysos.Utils.get_dims(_U_),
-        _X_,
-        _U_,
-    )
+    if model.time_type == CONTINUOUS
+        return MathematicalSystems.ConstrainedBlackBoxControlContinuousSystem(
+            dynamic(model, x_idx, u_idx),
+            Dionysos.Utils.get_dims(_X_),
+            Dionysos.Utils.get_dims(_U_),
+            _X_,
+            _U_,
+        )
+    else
+        return MathematicalSystems.ConstrainedBlackBoxControlDiscreteSystem(
+            dynamic(model, x_idx, u_idx),
+            Dionysos.Utils.get_dims(_X_),
+            Dionysos.Utils.get_dims(_U_),
+            _X_,
+            _U_,
+        )
+    end
 end
 
 function problem(model::Optimizer)
@@ -400,10 +429,13 @@ function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, value)
     return MOI.set(model.inner, attr, value)
 end
 
-export ∂, final, start
+export ∂, Δ, final, start
 
 function _diff end
-∂ = JuMP.NonlinearOperator(_diff, :diff)
+∂ = JuMP.NonlinearOperator(_diff, :∂)
+
+function _delta end
+Δ = JuMP.NonlinearOperator(_delta, :Δ)
 
 function _final end
 final = JuMP.NonlinearOperator(_final, :final)
@@ -411,6 +443,9 @@ final = JuMP.NonlinearOperator(_final, :final)
 function _start end
 start = JuMP.NonlinearOperator(_start, :start)
 
+function rem end
+rem_op = JuMP.NonlinearOperator(rem, :rem)
+
 # Type piracy
 function JuMP.parse_constraint_call(
     error_fn::Function,
@@ -429,3 +464,7 @@ function JuMP.parse_constraint_call(
     build_call = :(build_constraint($error_fn, $f, $(OuterSet)($set)))
     return parse_code, build_call
 end
+
+function Base.rem(x::JuMP.AbstractJuMPScalar, d)
+    return rem_op(x, d)
+end
diff --git a/src/optim/abstraction/uniform_grid_abstraction.jl b/src/optim/abstraction/uniform_grid_abstraction.jl
index f0021a335..37e53231c 100644
--- a/src/optim/abstraction/uniform_grid_abstraction.jl
+++ b/src/optim/abstraction/uniform_grid_abstraction.jl
@@ -4,6 +4,8 @@ module UniformGridAbstraction
 
 import StaticArrays
 
+import MathematicalSystems
+
 import Dionysos
 const DI = Dionysos
 const UT = DI.Utils
@@ -12,7 +14,7 @@ const ST = DI.System
 const SY = DI.Symbolic
 const PR = DI.Problem
 
-@enum ApproxMode GROWTH LINEARIZED DELTA_GAS DISCRETE_TIME
+@enum ApproxMode GROWTH LINEARIZED DELTA_GAS
 
 using JuMP
 
@@ -39,8 +41,7 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer
     approx_mode::ApproxMode
 
     ## for the discrete time (no need for the jacobian_bound) 
-    #but for the system_map, growthbound_map and sys_inv_map
-    system_map::Union{Nothing, Function}
+    # but for the `growthbound_map` and `sys_inv_map`
     growthbound_map::Union{Nothing, Function}
     sys_inv_map::Union{Nothing, Function}
 
@@ -63,7 +64,6 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer
             GROWTH,
             nothing,
             nothing,
-            nothing,
         )
     end
 end
@@ -97,23 +97,7 @@ function build_abstraction(concrete_system, model)
         error("Please set the `time_step`.")
     end
 
-    model.discretized_system = if model.approx_mode == GROWTH
-        if isnothing(model.jacobian_bound)
-            error("Please set the `jacobian_bound`.")
-        end
-        Dionysos.System.NewControlSystemGrowthRK4(
-            model.time_step,
-            concrete_system.f,
-            model.jacobian_bound,
-            noise,
-            noise,
-            model.num_sub_steps_system_map,
-            model.num_sub_steps_growth_bound,
-        )
-    elseif model.approx_mode == DISCRETE_TIME
-        if isnothing(model.system_map)
-            error("Please set the `system_map`.")
-        end
+    model.discretized_system = if concrete_system isa MathematicalSystems.ConstrainedBlackBoxControlDiscreteSystem
         if isnothing(model.growthbound_map)
             error("Please set the `growthbound_map`.")
         end
@@ -124,10 +108,23 @@ function build_abstraction(concrete_system, model)
             model.time_step,
             noise,
             noise,
-            model.system_map,
+            concrete_system.f,
             model.growthbound_map,
             model.sys_inv_map,
         )
+    elseif model.approx_mode == GROWTH
+        if isnothing(model.jacobian_bound)
+            error("Please set the `jacobian_bound`.")
+        end
+        Dionysos.System.NewControlSystemGrowthRK4(
+            model.time_step,
+            concrete_system.f,
+            model.jacobian_bound,
+            noise,
+            noise,
+            model.num_sub_steps_system_map,
+            model.num_sub_steps_growth_bound,
+        )
     elseif model.approx_mode == LINEARIZED
         Dionysos.System.NewControlSystemLinearizedRK4(
             model.time_step,
diff --git a/src/system/controlsystem.jl b/src/system/controlsystem.jl
index cd5acca73..ba23fdc58 100644
--- a/src/system/controlsystem.jl
+++ b/src/system/controlsystem.jl
@@ -276,7 +276,7 @@ end
 ############################################
 ############################################
 
-struct SymbolicSystem{}
+struct SymbolicSystem
     fsymbolicT::Any
     fsymbolic::Any
     Ts::Any