Control

BREAKING.md

@@ -4,6 +4,53 @@
 
 This document describes breaking changes in CTModels releases and how to migrate your code.
 
+## [0.9.8-beta] - 2026-03-16
+
+**No breaking changes** - This release adds piecewise constant interpolation for control signals while maintaining full backward compatibility.
+
+### New Features (Non-Breaking) - 0.9.8-beta
+
+- **Piecewise Constant Interpolation**
+  - New `ctinterpolate_constant` function with right-continuous steppost behavior
+  - Controls use `interpolation=:constant` by default in `build_solution`
+  - Control plotting uses `seriestype=:steppost` by default
+  - Enhanced `build_interpolated_function` with `interpolation` parameter
+
+- **Performance Improvements**
+  - Manual interpolation implementation ~20x-8600x faster to create
+  - 10-21% faster for multiple evaluations
+  - Zero allocations for interpolation object creation
+
+### API Enhancements (Non-Breaking)
+
+```julia
+# New constant interpolation (optional)
+interp = CTModels.ctinterpolate_constant(x, f)
+
+# Enhanced interpolation helpers (backward compatible)
+fu = OCP.build_interpolated_function(U, T, dim, type; interpolation=:constant)
+
+# Control plotting improvements (automatic)
+plot(sol, :control)  # Now uses seriestype=:steppost by default
+```
+
+### Migration Notes - 0.9.8-beta
+
+**No action required** - existing code continues to work unchanged.
+
+You can now benefit from improved control interpolation:
+
+```julia
+# Existing code (still works)
+sol = build_solution(ocp, T, X, U, v, P; objective=obj, ...)
+
+# New behavior (automatic)
+u = control(sol)  # Now uses piecewise constant interpolation
+plot(sol, :control)  # Now uses steppost plotting by default
+```
+
+---
+
 ## [0.9.7] - 2026-03-11
 
 **No breaking changes** - This release adds support for optimal control problems without a control input (`control_dimension == 0`) while maintaining full backward compatibility.

CHANGELOG.md

@@ -7,6 +7,59 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.9.8-beta] - 2026-03-16
+
+### 🚀 Major Features
+
+#### Piecewise Constant Interpolation for Control Signals
+
+- **New `ctinterpolate_constant` function**: Implements right-continuous steppost behavior for control signals
+- **Control interpolation**: Controls now use `interpolation=:constant` by default in `build_solution`
+- **Plotting integration**: Default `seriestype=:steppost` for control plotting, consistent with interpolation
+- **Performance optimized**: Manual implementation ~20x-8600x faster to create, 10-21% faster for multiple evaluations
+
+#### Enhanced Interpolation System
+
+- **Parameterized interpolation**: `build_interpolated_function` now accepts `interpolation::Symbol` (`:linear`, `:constant`)
+- **Manual `ctinterpolate`**: Replaced Interpolations.jl dependency with high-performance manual implementation
+- **Flat extrapolation**: Both interpolation functions use flat extrapolation (returns boundary values)
+
+### 📊 Performance Improvements
+
+- **Creation speed**: Manual interpolation objects are ~20x-8600x faster to create
+- **Evaluation speed**: 10-21% faster for multiple evaluations
+- **Memory efficiency**: Zero allocations for interpolation object creation
+- **Benchmark verified**: Comprehensive performance testing in `.extras/benchmark_interpolation.jl`
+
+### 🧪 Testing
+
+- **3456 tests pass**: Complete test coverage including new interpolation functionality
+- **75 interpolation tests**: Unit tests + comprehensive integration tests
+- **Behavior verification**: Tests confirm right-continuous steppost behavior
+- **Integration testing**: End-to-end testing from `build_solution` to `control()` interpolation
+- **Performance benchmarking**: Comprehensive testing in `.extras/benchmark_interpolation.jl`
+
+### 📝 API Changes
+
+```julia
+# New constant interpolation
+interp = CTModels.ctinterpolate_constant(x, f)  # Right-continuous steppost
+
+# Enhanced interpolation helpers
+fu = OCP.build_interpolated_function(U, T, dim, type; interpolation=:constant)
+
+# Control plotting with steppost (automatic)
+plot(sol, :control)  # Uses seriestype=:steppost by default
+```
+
+### 🔧 Internal Changes
+
+- **Dependency reduction**: Manual interpolation implementation reduces Interpolations.jl dependency
+- **Code clarity**: Explicit interpolation behavior with comprehensive documentation
+- **Cohesion**: Interpolation and plotting behaviors are now fully consistent
+
+---
+
 ## [0.9.7] - 2026-03-11
 
 ### Added - 0.9.7

Project.toml

@@ -1,12 +1,11 @@
 name = "CTModels"
 uuid = "34c4fa32-2049-4079-8329-de33c2a22e2d"
-version = "0.9.7"
+version = "0.9.8-beta"
 authors = ["Olivier Cots <olivier.cots@toulouse-inp.fr>"]
 
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
@@ -24,26 +23,10 @@ CTModelsJLD = "JLD2"
 CTModelsJSON = "JSON3"
 CTModelsPlots = "Plots"
 
-[extras]
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = [
-    "Aqua",
-    "JLD2",
-    "JSON3",
-    "Plots",
-    "Random",
-    "Test"
-]
-
 [compat]
 Aqua = "0.8"
 CTBase = "0.18"
 DocStringExtensions = "0.9"
-Interpolations = "0.16"
 JLD2 = "0.6"
 JSON3 = "1"
 LinearAlgebra = "1"
@@ -56,3 +39,18 @@ Random = "1"
 RecipesBase = "1"
 Test = "1"
 julia = "1.10"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = [
+    "Aqua",
+    "JLD2",
+    "JSON3",
+    "Plots",
+    "Random",
+    "Test"
+]

ext/plot.jl

@@ -96,8 +96,16 @@ function __plot_time!(
 
     #
     f(; kwargs...) = kwargs
+
+    # Default seriestype for controls (user can override with kwargs)
+    default_seriestype = if s == :control || s == :control_norm
+        :steppost
+    else
+        :path
+    end
+
     kwargs_plot = if isnothing(color)
-        f(; ylims=:auto, xlabel=t_label, ylabel=y_label, linewidth=2, z_order=:front, kwargs...)
+        f(; ylims=:auto, xlabel=t_label, ylabel=y_label, linewidth=2, z_order=:front, seriestype=default_seriestype, kwargs...)
     else
         f(;
             color=color,
@@ -106,6 +114,7 @@ function __plot_time!(
             ylabel=y_label,
             linewidth=2,
             z_order=:front,
+            seriestype=default_seriestype,
             kwargs...,
         )
     end

src/OCP/Building/interpolation_helpers.jl

@@ -7,7 +7,8 @@
 
 """
     _interpolate_from_data(data, T, dim, type_param; allow_nothing=false, 
-                           constant_if_two_points=false, expected_dim=nothing)
+                           constant_if_two_points=false, expected_dim=nothing,
+                           interpolation=:linear)
 
 Internal helper to create an interpolated function from discrete data.
 
@@ -19,6 +20,7 @@ Internal helper to create an interpolated function from discrete data.
 - `allow_nothing`: If false, throws IncorrectArgument when data is nothing
 - `constant_if_two_points`: If true and length(T)==2, return constant function
 - `expected_dim`: If provided, validates matrix dimension matches (via @ensure)
+- `interpolation`: Interpolation type (`:linear` or `:constant`)
 
 # Returns
 - Interpolated function (or nothing if data=nothing and allow_nothing=true)
@@ -38,6 +40,7 @@ function _interpolate_from_data(
     allow_nothing::Bool=false,
     constant_if_two_points::Bool=false,
     expected_dim::Union{Int,Nothing}=nothing,
+    interpolation::Symbol=:linear,
 )
     # Validation: nothing handling
     if isnothing(data)
@@ -83,7 +86,23 @@ function _interpolate_from_data(
     N = size(data, 1)
     cols = isnothing(dim) ? (:) : (1:dim)
     V = matrix2vec(data[:, cols], 1)
-    return ctinterpolate(T[1:N], V)
+
+    # Choose interpolation method
+    if interpolation == :linear
+        return ctinterpolate(T[1:N], V)
+    elseif interpolation == :constant
+        return ctinterpolate_constant(T[1:N], V)
+    else
+        throw(
+            Exceptions.IncorrectArgument(
+                "Invalid interpolation type";
+                got="interpolation=$interpolation",
+                expected=":linear or :constant",
+                suggestion="Use interpolation=:linear for linear interpolation or interpolation=:constant for piecewise-constant",
+                context="_interpolate_from_data",
+            ),
+        )
+    end
 end
 
 """
@@ -126,7 +145,8 @@ end
 
 """
     build_interpolated_function(data, T, dim, type_param; allow_nothing=false,
-                                constant_if_two_points=false, expected_dim=nothing)
+                                constant_if_two_points=false, expected_dim=nothing,
+                                interpolation=:linear)
 
 Unified function to build an interpolated function with deepcopy and scalar wrapping.
 
@@ -140,6 +160,7 @@ This is the main entry point that combines interpolation and wrapping in one cal
 - `allow_nothing`: Allow data=nothing (for optional duals)
 - `constant_if_two_points`: Return constant function if length(T)==2 (for costate)
 - `expected_dim`: Validate matrix has this dimension (for robustness)
+- `interpolation`: Interpolation type (`:linear` or `:constant`)
 
 # Returns
 - Wrapped interpolated function ready for use in Solution
@@ -154,6 +175,9 @@ This is the main entry point that combines interpolation and wrapping in one cal
 # State interpolation (required, with validation)
 fx = build_interpolated_function(X, T, dim_x, TX; expected_dim=dim_x)
 
+# Control with piecewise-constant interpolation
+fu = build_interpolated_function(U, T, dim_u, TU; expected_dim=dim_u, interpolation=:constant)
+
 # Costate with special 2-point handling
 fp = build_interpolated_function(P, T, dim_x, TP; 
                                  constant_if_two_points=true, expected_dim=dim_x)
@@ -172,6 +196,7 @@ function build_interpolated_function(
     allow_nothing::Bool=false,
     constant_if_two_points::Bool=false,
     expected_dim::Union{Int,Nothing}=nothing,
+    interpolation::Symbol=:linear,
 )
     # Handle T=nothing case
     if isnothing(T)
@@ -200,6 +225,7 @@ function build_interpolated_function(
         allow_nothing=allow_nothing,
         constant_if_two_points=constant_if_two_points,
         expected_dim=expected_dim,
+        interpolation=interpolation,
     )
 
     # Step 2: Wrap with deepcopy and scalar extraction

src/OCP/Building/solution.jl

@@ -235,8 +235,9 @@ function build_solution(
     # Build interpolated functions for state, control, and costate
     # Using unified API with validation and deepcopy+scalar wrapping
     # Note: costate uses its own grid (T_costate)
+    # Note: control uses piecewise-constant interpolation (steppost behavior)
     fx = build_interpolated_function(X, T_state, dim_x, TX; expected_dim=dim_x)
-    fu = build_interpolated_function(U, T_control, dim_u, TU; expected_dim=dim_u)
+    fu = build_interpolated_function(U, T_control, dim_u, TU; expected_dim=dim_u, interpolation=:constant)
     fp = build_interpolated_function(
         P, T_costate, dim_x, TP; constant_if_two_points=true, expected_dim=dim_x
     )
@@ -272,19 +273,22 @@ function build_solution(
         allow_nothing=true,
     )
     # Control box constraint duals share the control grid (T_control)
+    # Note: use piecewise-constant interpolation like control (steppost behavior)
     fccbd = build_interpolated_function(
         control_constraints_lb_dual,
         T_control,
         dim_control_constraints_box(ocp),
         Union{Matrix{Float64},Nothing};
         allow_nothing=true,
+        interpolation=:constant,
     )
     fccud = build_interpolated_function(
         control_constraints_ub_dual,
         T_control,
         dim_control_constraints_box(ocp),
         Union{Matrix{Float64},Nothing};
         allow_nothing=true,
+        interpolation=:constant,
     )
 
     # build Models

src/OCP/OCP.jl

@@ -43,8 +43,8 @@ include("aliases.jl")
 # Import macro from Utils module
 import ..Utils: @ensure
 
-# Import matrix2vec, ctinterpolate and to_out_of_place from Utils for solution building
-import ..Utils: matrix2vec, ctinterpolate, to_out_of_place
+# Import matrix2vec, ctinterpolate, ctinterpolate_constant and to_out_of_place from Utils for solution building
+import ..Utils: matrix2vec, ctinterpolate, ctinterpolate_constant, to_out_of_place
 
 # Load types first (no dependencies)
 include("Types/components.jl")

src/Utils/Utils.jl

@@ -11,6 +11,7 @@ including interpolation, matrix operations, and function transformations.
 The following functions are exported and accessible as `CTModels.function_name()`:
 
 - `ctinterpolate`: Linear interpolation for data
+- `ctinterpolate_constant`: Piecewise-constant interpolation for data
 - `matrix2vec`: Convert matrices to vectors
 
 # Private API
@@ -25,7 +26,6 @@ See also: `CTModels`
 module Utils
 
 using DocStringExtensions
-using Interpolations
 using CTBase: ctNumber
 
 # Private utilities (not exported)
@@ -37,6 +37,6 @@ include("interpolation.jl")
 include("matrix_utils.jl")
 
 # Export public API
-export ctinterpolate, matrix2vec
+export ctinterpolate, ctinterpolate_constant, matrix2vec
 
 end