SOLPS2ctrl.jl

find_files_in_allowed_folders(
    input_dirs::String...;
    eqdsk_file::String,
    recursive::Bool=true,
)

Searches a list of allowed folders for a set of filenames that will provide information about the SOLPS case. Returns a list of filenames with complete paths.

Example:

SOLPS2ctrl.find_files_in_allowed_folders(
    "<your samples folder>/D3D_Ma_184833_03600";
    eqdsk_file="g184833.03600",
)

geqdsk_to_imas!(
    eqdsk_file::String,
    dd::IMAS.dd;
    set_time::Union{Nothing, Float64}=nothing,
    time_index::Int=1,
)

Transfers the equilibrium reconstruction from an EFIT-style gEQDSK file into the IMAS DD structure.

preparation(
    eqdsk_file::String,
    dirs::String...;
    core_method::String="simple",
    filename::String="sd_input_data",
    output_format::String="json",
    eqdsk_set_time::Union{Nothing, Float64}=nothing,
    eq_time_index::Int=1,
    allow_boundary_flux_correction::Bool=false,
)::IMAS.dd

Gathers SOLPS and EFIT files and loads them into IMAS structure. Extrapolates profiles as needed to get a complete picture.

function add_rho_to_equilibrium(dd:IMAS.dd)

Adds equilibrium rho profile to the DD

check_rho_1d(
    dd::IMAS.dd;
    time_slice::Int=1,
    throw_on_fail::Bool=false,
)::Bool

Checks to see if rho exists and is valid in the equilibrium 1d profiles

extrapolate_core(
    edge_rho::Vector{Float64},
    edge_quantity::Vector{Float64},
    rho_output::Vector{Float64},
)::Vector{Float64}

Function for assuming a core profile when given edge profile data.

Concept:

1. value and derivative should be continuous when joining real data
2. derivative at magnetic axis is known to be 0 when making profile vs. rho, by the definition of rho
3. derivative probably does something fancier between the pedestal and axis than just linear interpolation, so add an extra point in there
4. there's a joint between the steep pedestal and the shallow core that needs an extra knot to manage it properly
5. after making up a very simple gradient profile out of a few line segments, integrate it to get the profile of the quantity in question

fill_in_extrapolated_core_profile!(
    dd::IMAS.dd,
    quantity_name::String;
    method::String="simple",
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    value_field::Symbol=:values,
    grid_ggd_idx::Int=1,
    cell_subset_idx::Int=5,
)

This function accepts a DD that should be populated with equilibrium and edge_profiles as well as a request for a quantity to extrapolate into the core. It then maps edge_profiles data to rho, calls the function that performs the extrapolation (which is not a simple linear extrapolation but has some trickery to attempt to make a somewhat convincing profile shape), and writes the result to core_profiles. This involves a bunch of interpolations and stuff.

Input arguments:

• dd: an IMAS data dictionary
• quantity_name: the name of a quantity in edge_profiles.profiles_2d and core_profiles.profiles_1d, such as "electrons.density"
• method: Extrapolation method.
• eq_time_idx: index of the equilibrium time slice to use. For a typical SOLPS run, the SOLPS mesh will be based on the equilibrium reconstruction at a single time, so the DD associated with the SOLPS run only needs one equilibrium time slice to be loaded. However, one could combine the complete equilibrium time series with the SOLPS run and then have to specify which slice of the equilibrium corresponds to the SOLPS mesh.
• eq_profiles_2d_idx: index of the profiles_2D in equilibrium time_slice.
• value_field: Symbolic name of the values field in quantity.
• grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS grid is fixed, so this index defaults to 1. But in future, if a time varying grid is used, then this index will need to be specified.
• cell_subset_idx: index of the subset of cells to use for the extrapolation. The default is 5, which is the subset of all cells. If edge_profiles data is instead present for a different subset, for instance, -5, which are b2.5 cells only, then this index should be set to -5.

mesh_psi_spacing(
    dd::IMAS.dd;
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    grid_ggd_idx::Int=1,
    space_idx::Int=1,
    avoid_guard_cell::Bool=true,
    spacing_rule="mean",
)

Inspects the mesh to see how far apart faces are in psi_N. Requires that GGD and equilibrium are populated.

Input Arguments:

• dd: a data dictionary instance with required data loaded into it
• eq_time_idx: index of the equilibrium time slice to use. For a typical SOLPS run, the SOLPS mesh will be based on the equilibrium reconstruction at a single time, so the DD associated with the SOLPS run only needs one equilibrium time slice to be loaded. However, one could combine the complete equilibrium time series with the SOLPS run and then have to specify which slice of the equilibrium corresponds to the SOLPS mesh.
• eq_profiles_2d_idx: index of the profiles_2D in equilibrium time_slice.
• grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS grid is fixed, so this index defaults to 1. But in future, if a time varying grid is used, then this index will need to be specified.
• space_idx: index of the space to use. For a typical SOLPS run, there will be only one space so this index will mostly remain at 1.
• avoid_guard_cell: assume that the last cell is a guard cell so take end-2 and end-1 instead of end and end-1
• spacing_rule: "edge" or "mean" to make spacing of new cells (in psi_N) be the same as the spacing at the edge of the mesh, or the same as the average spacing

cached_mesh_extension!(
    dd::IMAS.dd,
    eqdsk_file::String,
    b2fgmtry::String;
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    grid_ggd_idx::Int=1,
    space_idx::Int=1,
    clear_cache::Bool=false,
)::String

Adds an extended mesh to a data dictionary, possibly from a cached result.

Input Arguments:

• dd: The data dictionary. It will be modified in place.
• eqdsk_file: the name of the EQDSK file that was used to get equilibrium data in the dd.
• b2fgmtry: the name of the SOLPS geometry file that was used to get GGD info in edge_profiles in the dd.
• eq_time_idx: Index of the time slice in equilibrium
• eq_profiles_2d_idx: Index of the 2D profile set in equilibrium (there is usually only one)
• grid_ggd_idx: Index of the grid_ggd set in edge_profiles
• space_idx: Index of the space
• clear_cache: delete any existing cache file (for use in testing)

gas_unit_converter(
    value_in::Float64,
    units_in::String,
    units_out::String;
    species::String="H",
    temperature::Float64=293.15,
)

Converts gas flows between different units. Uses ideal gas law to convert between Pressure * volume type flows / quantities and count / current types of units. There is a version that accepts floats in and outputs floats, and another that deals in Unitful quantities.

gas_unit_converter(
    value_in::Unitful.Quantity,
    units_in::String,
    units_out::String;
    species::String="H",
    temperature=293.15 * Unitful.K,
)

Converts gas flows between different units. Uses ideal gas law to convert between Pressure * volume type flows / quantities and count / current types of units. This is the Unitful version.

Output will be unitful, but the units are not simplified automatically. You can perform operations such as

• (output |> Unitful.upreferred).val
• Unitful.uconvert(Unitful.whatever, output).val

to handle simplification or conversion of units.

Although this function pretends torr L s$^{-1}$ and Pa m$^3$ s$^{-1}$ are different, use of Unitful should cause them to behave the same way as long as you simplify or convert units at the end. This means that you can use other pressure*volume type gas units and call them torr L s$^{-1}$ and the script will deal with them up to having messy units in the output.

Plant

Abstract parent type for all plants. Whenever a user defined plant is created, it must be subtype of Plant.

To create a new plant with custom feautres, it must be defined as a mutable stucture which is daughter of Plant that contains all settings and state information for the plant and the instance itself should be callable to take as input a Vector{Float64} and should return an output of Vector{Float64}.

mutable struct CustomPlant <: Plant
    settings
    state
    # ... Anything more
end

function (cp::CustomPlant)(inp::Vector{Float64})::Vector{Float64}
    # perform the plant single step forward calcualtion with inp
    # update cp.state if required
    return output
end

NOTE: If you need to add a linear system in your plant model, add a LinearPlant instance in the attributes of your CustomPlant and just call the LinearPlant inside your function call.

LinearPlant

Implementation of linear system with option of scaling and offseting inputs and outputs. It stores the system in sys and state of the system in x. Constructor:

LinearPlant(
    sys::Union{PredictionStateSpace, StateSpace},
    x=zeros(Float64, size(sys.A, 1));
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
)

Creates a LinearPlant instance with state space system sys and state vector x. It also defined input offsetting and scaling with inp_offset and inp_factor, and similarly output unscaling and unoffseting with out_offset and out_factor.

InputConditioning

Abstract parent type for creating input conditioning to plants. Whenever a user defined input coniditioning is created it must be subtype of InputConditioning.

To create a new input conditioning with custom feautres, it must be defined as a mutable stucture which is daughter of InputConditioning that contains all settings and state information for the function and the instance itself should be callable to take as input a Vector{Float64} and kwargs for any parameters used and should return an output of Vector{Float64} which is same size as the input.

mutable struct CustomInpCond <: InputConditioning
    settings
    state
    # ... Anything more
end

function (inp_cond::CustomInpCond)(inp::Vector{Float64}; kwargs...)::Vector{Float64}
    # perform the input coniditioning single step forward calcualtion with inp
    # Use parameters from kwargs
    # update inp_cond.state if required
    return output
end

Note that inp_cond must have a call signature of inp_cond(inp::Vector{Float64}; kwargs...)::Vector{Float64} to take input as a vector for multiple inputs at a particular time instance

InpCondLinPlant

Implementation of a LinearPlant with an added input conditioning which could be used to make it non-linear. It stores the LinearPlant in linear_plant, the input coniditioning function in inp_cond and any keyword arguments for the inp_cond in inp_cond_kwargs. Constructor:

InpCondLinPlant{T}(
    linear_plant::LinearPlant,
    inp_cond::InputConditioning
    inp_cond_kwargs::Dict{Symbol, T}
) where {T}

Creates a InpCondLinPlant instance. inp_cond_kwargs can be used to have changeable parameters in the input conditioning which can be used for optimization and fiting.

system_id(
    inp::Union{Vector{Float64}, Matrix{Float64}},
    out::Union{Vector{Float64}, Matrix{Float64}},
    tt::Vector{Float64},
    order::Int;
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
    inp_cond::Union{Nothing, InputConditioning}=nothing,
    inp_cond_kwargs::Dict{Symbol, T}=Dict{Symbol, Any}(),
    newpem_kwargs::Dict{Symbol, U}=Dict{Symbol, Any}(),
    verbose::Bool=false,
) where {T, U}

Perform system identification for a set on input data inp, output data out, and time series vector tt. If there are more than one inputs or outputs, provide them as Matrix with first dimension for ports (input or output) and second dimension for time.

If inp_cond is provided, it is applied before offsetting and scaling for performing system identification with keywords for this function provided in inp_cond_kwargs.

This function uses ControlSystemIdentification.newpem to perform the system identification. Any additional keywords for this function should be passed as dictionary in newpem_kwargs. For advanced use, it is recommended to do system identification directly with newpem instead of using this function.

Returns a linear state space model of the order.

system_id_optimal_inp_cond(
    inp::Union{Vector{Float64}, Matrix{Float64}},
    out::Union{Vector{Float64}, Matrix{Float64}},
    tt::Vector{Float64},
    order::Int,
    inp_cond::InputConditioning,
    inp_cond_args_guess::Dict{Symbol, T};
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
    inp_cond_args_lower::Dict{Symbol, V}=Dict{Symbol, Any}(),
    inp_cond_args_upper::Dict{Symbol, W}=Dict{Symbol, Any}(),
    newpem_kwargs::Dict{Symbol, U}=Dict{Symbol, Any}(),
    curve_fit_kwargs::Dict{Symbol, X}=Dict{Symbol, Any}(),
    verbose::Bool=false,
) where {T, U, V, W}

Perform system identification for a set on input data inp, output data out, and time series vector tt. If there are more than one inputs or outputs, provide them as Matrix with first dimension for ports (input or output) and second dimension for time.

The inp_cond is applied before offsetting and scaling for performing system identification. The inp_cond_args_guess is used as initial keyword arguments that provide the parameters of the inp_cond. These arguments are then used to find the best fit while iteratively performing system_id in each step.

This function uses ControlSystemIdentification.newpem to perform the system identification. Any additional keywords for this function should be passed as dictionary in newpem_kwargs.

This function uses LsqFit.curve_fit to fit the parameters of input conditions along with performing the system identification.

For advanced use, it is recommended to do system identification directly with newpem and optimize using your favorite fitting method instead of using this function.

Returns a NonLinInpLinearPlant that has stored best fit parameters for the input non-linearity and a best fit linear model.

model_evolve(
    plant::Plant,
    inp::Union{Vector{Float64}, Matrix{Float64}};
    x0::Union{Nothing, Vector{Float64}}=nothing,
    initialize_x::Bool=false,
)

Evolve a plant model through a time series of inputs. If inp is a Matrix, the first dimension will hold different inputs and second dimension will be along time. If inp is a vector, it would be assumed that the model has a single input. If the model also happens to have a single output and inp is a vector, the returned outputs will be a vector as well. If x0 is not provided, the stored state vector of the plant model will be used for initialization. If initialize_x is true, then the state vector of the plant model would be initialized to match first input value so that no sudden jumps happen.

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}.

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 instance in the attributes of your CustomActuator and just call the DelayedActuator inside your function call.

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.

Controller

Abstract parent type for all controllers. Whenever a user defined controller is created it must be a subtype of Controller.

To create a new controller algorithm, it should be defined as a mutable structure that is a daughter of Controller and should contain all settings and state information to be stored. It must be a callable structure that can use any of the following keyword arguments:

• ii::Int: Iteration index
• target::Matrix{Float64}: Target waveform (No. of signals x Time steps)
• plant_inp::Matrix{Float64}: Inputs to plant (No. of inputs x Time steps)
• plant_out::Matrix{Float64}: Outputs from plant (No. of outputs x Time steps)
• act::Actuator: Actuator model in the loop
• inp_feedforward::Matrix{Float64}: Feedforward input to the plant (No. of inputs x Time steps)
• kwargs..: Required to ignore unused keyword arguments

LinearController

Implementation for using any linear controller. It stores ctrl_ss::StateSpace{TE} where {TE <: Discrete} for storing any linear controller as a discrete state space model using ControlSystemsBase.ss. It also stores the state vector for the state space model as ctrl_x0::Vector{Float64}. It's call signature is:

(lc::LinearController)(;
    ii::Int,
    target::Matrix{Float64},
    plant_out::Matrix{Float64},
    kwargs...,
)::Vector{Float64}

Calcualtes error as target[:, ii] .- plant_out[:, ii] and runs it through lsim_step().

PVLC

# Constructor
PVLC(
    ctrl_ss::StateSpace{TE},
    ctrl_x0::Vector{Float64},
    plant::LinearPlant,
    h::Int,
) where {TE <: Discrete}

Implementation of Predicted Variable Linear Controller (PVLC). It stores ctrl_ss::StateSpace{TE} where {TE <: Discrete} for storing any linear controller as a discrete state space model using ControlSystemsBase.ss. It also stores the state vector for the state space model as ctrl_x0::Vector{Float64}. Additionally, it stores the plant, the number of steps of history h used for state tracking and state prediction matrices Y2x and U2x from state_prediction_matrices().

This controller has a call signature:

(pvlc::PVLC)(;
    ii::Int,
    target::Matrix{Float64},
    plant_inp::Matrix{Float64},
    plant_out::Matrix{Float64},
    act::Actuator,
    kwargs...,
)::Vector{Float64}

Tracks the state vector of the plant using h steps of history from plant_inp and plant_out and uses it to calculate future output of the plant. It compares it with a future target value and applies the linear controller ctrl_ss there.

state_prediction_matrices(
    sys::Union{PredictionStateSpace, StateSpace},
    h::Int=size(sys.A, 1),
)

Calculate state prediction matrices for a linear system sys given h steps of input and output history. This function returns two matrices, Y2x and U2x that can be used to calculate least square fitted estimate of the current state vector of the system. Y2x has size: (State size of sys) x (No. of outputs times h) U2x has size: (State size of sys) x (No. of inputs times h) The estimated state vector is obtained by: Y2x * Y + U2x * U where Y is all the h outputs of system stacked into a single vector. and U is all the h inputs to the system stacked into a single vector.

MPC

# Constructor
MPC(
    plant::Plant,
    h::Int,
    act::Actuator,               # Actuator model without delay
    horizon::Int,                # Number of steps in future after latency
    nopt::Int,                   # Number of optimization points in horizon window
    opt_every::Int;              # Run cost optimization every opt_every steps
    ctrl_out_bounds::Tuple{Vector{Float64}, Vector{Float64}}=(
        Array{Float64}(undef, 0),
        Array{Float64}(undef, 0),
    ),
    guess::Union{Symbol, Vector{Float64}}=:zeros,
    curve_fit_kwargs::Dict{Symbol, T}=Dict{Symbol, Any}(),
)

Implementation of simple leaqt square optimized Model Predictive Controller (MPC). It stores the plant, the number of steps of history h used for state tracking and state prediction matrices Y2x and U2x from state_prediction_matrices(). It stores a current deep copy of actuator instance in act to try it during optimization and it stores setting for least square optimization. It also stores a future_evolve that is created based on plant, act, and optimization setting and is the model function that is used for optimization later. This controller stores a buffer for control outputs, ctrl_out_buffer so that it can be called less often and it can reuse it's previous optimization results.

This contructor takes in minimum required information to create a self-consistent MPC instance. It sets the other dependent quantities in MPC such as Y2x, U2x, min_delay, and create a future_evolve function and initializes the ctrl_out_buffer. act needs to be a deepcopy of the actuator instance. horizon is the number of steps after the min_delay among all acturators for which the optimization is carried out. nopt is the number of optimization points taken along the horizon which are lineary distributed. The gaps between this optimzation points are interpolated linearly in the output. opt_every defines the frequency of optimization, i.e. at every opt_every call of this controller, the optimization is carried out. This avoids unnecessary over-calculations and thus results in a faster controller. ctrl_out_bounds is a tuple of lower bounds and upper bounds for the control output. guess is used to create the initial guess during least square optimization. If :last, it would use the last controller output as initial setting. If it is a Vector, each initialization starts with this Vector. In all other cases, initial guess is zeros. curve_fit_kwargs can be used to provide keyword arguments that go to curve_fit.

This controller has a call signature:

function (mpc::MPC)(;
    ii::Int,
    target::Matrix{Float64},
    plant_inp::Matrix{Float64},
    plant_out::Matrix{Float64},
    act::Actuator,
    inp_feedforward::Matrix{Float64}=zeros(
        Float64,
        (size(get_sys(plant).B, 2), length(target)),
    ),
    kwargs...,
)::Vector{Float64} where {T}

Tracks the state vector of the plant using h steps of history from plant_inp and plant_out and uses it along with act to run an optimization to match target in future after the minimum delay from all the actuators. This function uses curve_fit() for the non-linear least squares fitting which uses Levenberg-Marquardt algorithm. This function performs the optimization every opt_every call and uses the stored control output in ctrl_out_buffer meanwhile.

run_closed_loop_sim(
    plant::Union{PredictionStateSpace{TE}, StateSpace{TE}},
    act::Actuator,
    ctrl::Controller,
    target::Matrix{Float64};
    inp_feedforward::Matrix{Float64}=zeros(
        Float64,
        (size(get_sys(plant).B, 2), size(target, 2)),
    ),
    ctrl_start_ind::Int=1,
    noise_plant_inp::Matrix{Float64}=zeros(
        Float64,
        (size(get_sys(plant).B, 2), size(target, 2)),
    ),
    noise_plant_out::Matrix{Float64}=zeros(
        Float64,
        (size(get_sys(plant).C, 1), size(target, 2)),
    ),
    noise_ctrl_out::Matrix{Float64}=zeros(
        Float64,
        (size(get_sys(plant).B, 2), size(target, 2)),
    ),
) where {T, TE <: Discrete}

Generic function to run closed loop simulations with provided plant, actuator act, controller ctrl, and target waveform target. The length of simulation is determined by provided target. Keyword arguments are possible for providing adjustments to inputs and outputs of the plant model as explained in model_evolve(). Additionally, ctrl_start_ind can be provided to start control loop at an arbitrary point in the loop. noise_plant_inp, noise_plant_out, and noise_ctrl_out allow addition of predefined noise waveforms at the input of plant, output of plant, and the output of controller respectively.

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

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

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