Reorganization and introducing Plant as a type

Project.toml

@@ -1,7 +1,7 @@
 name = "SOLPS2ctrl"
 uuid = "a531d12f-ac8a-43e8-b6d9-bd121431dd49"
 authors = ["David Eldon <eldond@fusion.gat.com>"]
-version = "2.2.1"
+version = "2.2.2"
 
 [deps]
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"

docs/src/index.md

@@ -30,15 +30,16 @@ using Pkg
 Pkg.add("SOLPS2ctrl")
 ```
 
-## Top file handling functions
+## SOLPS related functionality
+### Top file handling functions
 
 ```@docs
 find_files_in_allowed_folders
 geqdsk_to_imas!
 preparation
 ```
 
-## Repairing/filling out partial equilibrium files
+### Repairing/filling out partial equilibrium files
 
 Tools for repairing/filling out partial equilibrium files.
 
@@ -51,18 +52,18 @@ add_rho_to_equilibrium!
 check_rho_1d
 ```
 
-## Extrapolations
+### Extrapolations
 
 Utilities for extrapolating profiles
 
-### Core profile extrapolations
+#### Core profile extrapolations
 
 ```@docs
 extrapolate_core
 fill_in_extrapolated_core_profile!
 ```
 
-### Edge profiles extrapolations
+#### Edge profiles extrapolations
 
 These functions have not been fully tested and/or supported yet.
 
@@ -71,17 +72,38 @@ mesh_psi_spacing
 cached_mesh_extension!
 ```
 
-## System identification and modeling
+### Unit conversion utilities
+
+```@docs
+gas_unit_converter
+```
+
+## Control related functionality
+
+This package provides a platform to perform closed loop simulations of controllers with
+plant models and actuator models. The key function to run these simulations is `run_closed_loop_sim()` described below. The common architecture for this platform is to create `mutable struct` for each of the actuator, plant, and controller which are
+subtypes of provided abstract types `Actuator`, `Plant`, and `Controller` and the
+struct itself is callable and performs the required function for a discrete time
+instance. See test examples to get an idea of how to use this feature.
+
+### Plant
+
+```@docs
+Plant
+LinearPlant
+InputConditioning
+InpCondLinPlant
+```
+
+### System identification and modeling
 
 ```@docs
-offset_scale
-unscale_unoffset
 system_id
 system_id_optimal_inp_cond
 model_evolve
 ```
 
-## Actuators
+### Actuators
 
 ```@docs
 Actuator
@@ -98,26 +120,26 @@ for (i, inp) in enumerate(inp_series)
 end
 ```
 
-## Controllers
+### Controllers
 
 ```@docs
 Controller
 ```
 
-### Linear Controller
+#### Linear Controller
 
 ```@docs
 LinearController
 ```
 
-### Predicted Variable Linear Controller
+#### Predicted Variable Linear Controller
 
 ```@docs
 PVLC
 state_prediction_matrices
 ```
 
-#### State Prediction Matrix Algebra
+##### State Prediction Matrix Algebra
 
 At step k:
 
@@ -250,22 +272,22 @@ Then, future state can be predicted by:
 x_{k+1} = \mathcal{N} \mathcal{L}^{-1} (\vec{Y} - \mathcal{M} \vec{U}) + \mathcal{O} \vec{U}
 ```
 
-### Model Predictive Controller
+#### Model Predictive Controller
 
 ```@docs
 MPC
 ```
 
-## Closed loop simulations
+### Closed loop simulations
 
 ```@docs
-lsim_step
-model_step
 run_closed_loop_sim
 ```
 
-## Unit conversion utilities
+### Control related utilities
 
 ```@docs
-gas_unit_converter
-```
+lsim_step
+offset_scale
+unscale_unoffset
+```

src/SOLPS2ctrl.jl

@@ -10,7 +10,10 @@ export find_files_in_allowed_folders, geqdsk_to_imas!, preparation
 include("$(@__DIR__)/supersize_profile.jl")
 include("$(@__DIR__)/repair_eq.jl")
 include("$(@__DIR__)/unit_utils.jl")
+include("$(@__DIR__)/control_utils.jl")
+include("$(@__DIR__)/plant.jl")
 include("$(@__DIR__)/system_id.jl")
+include("$(@__DIR__)/actuators.jl")
 include("$(@__DIR__)/controllers.jl")
 
 """

src/actuators.jl

@@ -0,0 +1,100 @@
+import DataStructures: Queue, enqueue!, dequeue!
+
+export Actuator, DelayedActuator
+
+"""
+    Actuator
+
+Abstract parent type for all actuators. Whenever a user defined actuator is created,
+it must be subtype of `Actuator`.
+
+To create a new actuator with custom function, it must be defined as a mutable stucture
+which is daughter of Actuator that contains all settings and state information for
+the actuator and the instance itself should be callable to take as input a
+`Vector{Float64}` and should return an output of `Vector{Float64}`.
+
+```julia
+mutable struct CustomActuator <: Actuator
+    settings
+    state
+    # ... Anything more
+end
+
+function (ca::CustomActuator)(inp::Vector{Float64})::Vector{Float64}
+    # perform the actuation calcualtions with inp
+    # update ca.state if required
+    return output
+end
+```
+
+**NOTE:** If you need to add a delay in the actuator, add a [`DelayedActuator`](@ref)
+instance in the attributes of your `CustomActuator` and just call the DelayedActuator
+inside your function call.
+"""
+abstract type Actuator end
+
+"""
+    DelayedActuator{U}
+
+Implementation of delayed actuation. It stores `delay::Int` for number of time steps of
+delay and `buffer::Queue{U}` which stores the delayed actuations in a queue.
+Constructor:
+
+    DelayedActuator(
+        delay::Int;
+        default_output::T=[0.0],
+    ) where {T}
+
+Creates a DelayedActuator{T} instance with `delay` and initializes the `buffer`
+pre-filled upto brim with the default_output.
+"""
+mutable struct DelayedActuator{U} <: Actuator
+    delay::Int
+    buffer::Queue{U}
+    function DelayedActuator(
+        delay::Int;
+        default_output::T=[0.0],
+    ) where {T}
+        buffer = Queue{T}(delay)
+        for i ∈ 1:delay
+            enqueue!(buffer, default_output)
+        end
+        return new{T}(delay, buffer)
+    end
+end
+
+function (da::DelayedActuator)(inp::Union{Float64, Vector{Float64}})::Vector{Float64}
+    if da.delay > 0
+        enqueue!(da.buffer, inp)
+        return dequeue!(da.buffer)
+    else
+        return inp
+    end
+end
+
+function get_future_inp(act::Actuator)
+    for fieldname ∈ fieldnames(typeof(act))
+        if fieldtype(typeof(act), fieldname) == DelayedActuator
+            return get_future_inp(getfield(act, fieldname))
+        end
+    end
+    return Matrix{Float64}(undef, 0, 0)
+end
+
+get_future_inp(act::DelayedActuator) = stack(act.buffer)
+
+function get_min_delay(act::Actuator)::Int
+    min_delay = typemax(Int)
+    has_delay = false
+    for fieldname ∈ fieldnames(typeof(act))
+        if fieldtype(typeof(act), fieldname) == DelayedActuator
+            min_delay = min(min_delay, getfield(act, fieldname).delay)
+            has_delay = true
+        end
+    end
+    if has_delay
+        return min_delay
+    else
+        return 0
+    end
+end

src/control_utils.jl

@@ -0,0 +1,60 @@
+import ControlSystemIdentification: PredictionStateSpace
+import ControlSystemsBase: StateSpace
+
+export lsim_step, offset_scale, unscale_unoffset
+
+"""
+    lsim_step(
+        sys::Union{PredictionStateSpace, StateSpace}, u::Vector{Float64};
+        x0::Vector{Float64}=zeros(size(sys.A, 1)),
+    )::Tuple{Vector{Float64}, Vector{Float64}}
+
+Single step version of [ControlSystemsBase.lsim](https://juliacontrol.github.io/ControlSystems.jl/dev/lib/timefreqresponse/#ControlSystemsBase.lsim-Tuple%7BAbstractStateSpace,%20AbstractVecOrMat,%20AbstractVector%7D)
+"""
+function lsim_step(
+    sys::Union{PredictionStateSpace, StateSpace}, u::Vector{Float64};
+    x0::Vector{Float64}=zeros(size(sys.A, 1)),
+)::Tuple{Vector{Float64}, Vector{Float64}}
+    x = sys.A * x0 + sys.B * u
+    y = sys.C * x0 + sys.D * u
+    return y, x
+end
+
+"""
+    offset_scale(
+        val::Union{Float64, Vector{Float64}, Matrix{Float64}};
+        offset::Union{Float64, Vector{Float64}}=0.0,
+        factor::Union{Float64, Vector{Float64}}=1.0,
+    )::typeof(val)
+
+Subtract an `offset` and multiply by a `factor`, the `val` to make it nominally in the
+range of -1 to 1 (not strictly) for easy identification of system.
+
+    (val .- offset) .* factor
+"""
+function offset_scale(
+    val::Union{Float64, Vector{Float64}, Matrix{Float64}};
+    offset::Union{Float64, Vector{Float64}}=0.0,
+    factor::Union{Float64, Vector{Float64}}=1.0,
+)::typeof(val)
+    return (val .- offset) .* factor
+end
+
+"""
+    unscale_unoffset(
+        offset_scaled::Union{Float64, Vector{Float64}, Matrix{Float64}};
+        offset::Union{Float64, Vector{Float64}}=0.0,
+        factor::Union{Float64, Vector{Float64}}=1.0,
+    )::typeof(offset_scaled)
+
+Undo previously applied offset and scaling.
+
+    offset_scaled ./ factor .+ offset
+"""
+function unscale_unoffset(
+    offset_scaled::Union{Float64, Vector{Float64}, Matrix{Float64}};
+    offset::Union{Float64, Vector{Float64}}=0.0,
+    factor::Union{Float64, Vector{Float64}}=1.0,
+)::typeof(offset_scaled)
+    return offset_scaled ./ factor .+ offset
+end