Last updated: 2021-08-03

Checks: 7 0

Knit directory: VSA_altitude_hold/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20210617)

The command set.seed(20210617) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 6a83009

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 6a83009. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    renv/library/
    Ignored:    renv/local/
    Ignored:    renv/staging/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/simulation_framework.Rmd) and HTML (docs/simulation_framework.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
html	6a83009	Ross Gayler	2021-08-03	WIP
html	dcdf309	Ross Gayler	2021-07-20	WIP
Rmd	6c3ace9	Ross Gayler	2021-07-20	WIP
html	6c3ace9	Ross Gayler	2021-07-20	WIP

This notebook documents the functions used to simulate the repeated updating of the state of the PID controller.

The general issues around the software interface between the VSA and multicopter components are discussed in Relation of VSA and multicopter simulation components.

The specific software design issues around the simulation of the PID controller are discussed in Structure of the VSA data flow graph implementation.

1 `run_simulation()`

The function to run the simulation is defined as:

# Function to run an arbitrary simulation
run_simulation <- function(
  input_df, # dataframe[n_steps, n_input_vars] - values of all input variables at all times
  init_state, # dataframe[1, n_state_vars] - initial values of state variables used by f_update()
  f_update # function(prev_state, input) - state update function, args are 1-row dataframes
) # value - dataframe[n_steps, n_input_vars + n_state_vars + 1]
  # One row per time step
  # One column for each input variable, state variable, and time (t)
{
  # Apply the state update to tghe input values
  state_df <- input_df %>% 
    base::split(seq(nrow(.))) %>% # split input into a list of 1-row data frames
    purrr::accumulate( # accumulate list of simulated states
      f_update,
      .init = init_state
    ) %>% 
    purrr::discard(.p = seq_along(.) == 1) %>% # discard first element (initial state)
    dplyr::bind_rows() %>%  # convert list of time step states to a data frame
    dplyr::bind_cols(input_df, .) %>% # add input variables
    dplyr::mutate(t = 1:nrow(input_df), .before = everything()) # add time variable
}

# The input variables data frame is split into a list of one-row data frames
# because accumulate() requires it. The input variables could have been created
# in that format, but that's a much less convenient for general manipulation
# and the list of rows format is only required for accumulate().
# That's why the reformatting occurs on the fly in run_simulation().

# The time and input variables are added to the data frame for convenience.

The arguments are:

input_df - The history of the input variables.
This is a data frame with one row for each time step (\(t = 1 \ldots n_{steps}\)) and one column for each input variable. For this project the input variables are numeric scalars, but in general they can be arbitrarily complex structures if the column is a list column.
init_state - The initial state of the state variables.
This is a a one-row date frame. There is one column for each state variable. It must contain those state variables where the update function depends on the previous value of that variable. Those state variables which are VSA vectors must be list columns.
f_update - The state update function.
This is a function of two arguments. The first argument is a one-row data frame containing the previous values of the state variables. The second argument is a one-row data frame containing the input variables for that time step. The value of the state update function is a one-row data frame containing the updated values of the state variables.

The value of the run_simulation() function is a data frame with columns for the input variables, state variables, and time (\(t\)). There is one row for each of the time steps (\(t = 1 \ldots n_{steps}\)).

2 Set up example

The following sections set up an example of using run_simulation(). This is a drastically scaled down example in that the number of time steps and VSA vector dimensions are greatly reduced so that the simulation result can be viewed in its entirety.

2.1 Set global parameters

Set the dimensionality of the VSA vectors and the number of time steps in the simulation.

v_dim <- 10 # vector dimension of state variables (in practice = 1e4)
n_step <- 5 # number of time steps in simulation (in practice = 500)

2.2 Input variables

Create a data frame to hold the input values. Normally these will be from a file exported from the multicopter simulation.

The input values are contained in a data frame.

Each of the input variables is a time series - one value per time step. The data frame contains as many input variables as you need.

For this example we are just generating random input data. The input data for this project are all numeric scalars, so it could have been stored in a matrix, but in general, the inputs might be mixed types or complex structures so we use a data frame for generality.

input_df <- tibble::tibble(
  i_x = rnorm(n = n_step),
  i_y = rnorm(n = n_step)
)

input_df

# A tibble: 5 × 2
     i_x     i_y
   <dbl>   <dbl>
1 -0.437 -0.634 
2  0.329  0.628 
3 -0.447 -0.752 
4 -0.183  0.0409
5  1.67  -2.64

2.3 Initial state

The initial state consists of whatever state variables the state update function refers to as previous state.

The previous state is represented as a one-row data frame with one column for each of the required state variables. The VSA vector variables need to be stored in list columns.

I expect that we will generally use a function to construct the initial state.

# Function to make the initial state
mk_init_state <- function(
  v_dim, # integer[1] - dimension of vector-valued state variables
  seed = NULL # integer[1] - random generation seed
) # value - dataframe[1, ] - initial state of system
{
  set.seed(seed)
  
  # Construct initial state
  # MUST contain all columns used as previous state by update_state()
  # Any other state columns are optional
  tibble::tibble_row(
    # scalar valued state variables
    s_x = 0, s_y = 0,
    # vector valued state variables (each must be wrapped in list())
    s_0 = sample(c(-1, 1), size = v_dim, replace = TRUE) %>% list(), 
    s_1 = sample(c(-1, 1), size = v_dim, replace = TRUE) %>% list()
  )
}

# Make the initial state
init_state <- mk_init_state(v_dim, seed = 42)

init_state

# A tibble: 1 × 4
    s_x   s_y s_0        s_1       
  <dbl> <dbl> <list>     <list>    
1     0     0 <dbl [10]> <dbl [10]>

str(init_state)

tibble [1 × 4] (S3: tbl_df/tbl/data.frame)
 $ s_x: num 0
 $ s_y: num 0
 $ s_0:List of 1
  ..$ : num [1:10] -1 -1 -1 -1 1 1 1 1 -1 1
 $ s_1:List of 1
  ..$ : num [1:10] -1 1 -1 1 -1 -1 1 1 1 1

2.4 State update function

The state update function takes the previous state and current input variables as arguments and returns the current state.

This specific update function is a gibberish example, but is generically the sort of thing I will be doing

The previous state is given as a one-row data frame. The VSA vector values must be stored as list columns.

The current input variables are given as a one-row data frame. For this project the input variables are numeric scalars which are represented as standard columns, but more generally the inputs may be more complex and require list columns.

The value returned by the state update function is a one-row data frame containing the updated values of all the state variables.

# Function to update state from previous state and current input
update_state <- function(
  prev_state, # dataframe[1, ] - previous state
  input # dataframe[1, ] - current input
) 
{
  # Update scalar valued state variables
  s_x = prev_state$s_x + input$i_x 
  s_y = prev_state$s_y + input$i_y 
  
  # Update vector valued state variables
  # List columns must be unlisted before use
  s_0 = unlist(prev_state$s_0) * s_x
  s_1 = unlist(prev_state$s_1) * s_y
  s_2 = s_0 * s_1
  s_3 = s_2 * s_0
  s_4 = s_2 * s_1
  s_5 = s_3 + s_4
  
  # Construct and return a 1-row dataframe containing the updated state
  tibble::tibble_row(
    # Scalar valued state variables
    s_x, s_y, 
    # Vector valued state variables
    # List columns must be wrapped by list()
    s_0 = list(s_0),
    s_1 = list(s_1),
    s_2 = list(s_2),
    s_3 = list(s_3),
    s_4 = list(s_4),
    s_5 = list(s_5)
  )
}

3 Run small simulation

Run the simulation based on the previous settings.

# Run the simulation, returning the history as a data frame
history_df <- input_df %>% run_simulation(init_state, update_state)

# Take a look at the shape of the history
history_df

# A tibble: 5 × 11
      t    i_x     i_y    s_x      s_y s_0        s_1    s_2   s_3   s_4   s_5  
  <int>  <dbl>   <dbl>  <dbl>    <dbl> <list>     <list> <lis> <lis> <lis> <lis>
1     1 -0.437 -0.634  -0.437 -0.634   <dbl [10]> <dbl … <dbl… <dbl… <dbl… <dbl…
2     2  0.329  0.628  -0.109 -0.00641 <dbl [10]> <dbl … <dbl… <dbl… <dbl… <dbl…
3     3 -0.447 -0.752  -0.555 -0.759   <dbl [10]> <dbl … <dbl… <dbl… <dbl… <dbl…
4     4 -0.183  0.0409 -0.738 -0.718   <dbl [10]> <dbl … <dbl… <dbl… <dbl… <dbl…
5     5  1.67  -2.64    0.933 -3.36    <dbl [10]> <dbl … <dbl… <dbl… <dbl… <dbl…

str(history_df)

tibble [5 × 11] (S3: tbl_df/tbl/data.frame)
 $ t  : int [1:5] 1 2 3 4 5
 $ i_x: num [1:5] -0.437 0.329 -0.447 -0.183 1.672
 $ i_y: num [1:5] -0.6343 0.6279 -0.7524 0.0409 -2.6394
 $ s_x: num [1:5] -0.437 -0.109 -0.555 -0.738 0.933
 $ s_y: num [1:5] -0.63429 -0.00641 -0.75879 -0.71792 -3.35736
 $ s_0:List of 5
  ..$ : num [1:10] 0.437 0.437 0.437 0.437 -0.437 ...
  ..$ : num [1:10] -0.0475 -0.0475 -0.0475 -0.0475 0.0475 ...
  ..$ : num [1:10] 0.0264 0.0264 0.0264 0.0264 -0.0264 ...
  ..$ : num [1:10] -0.0195 -0.0195 -0.0195 -0.0195 0.0195 ...
  ..$ : num [1:10] -0.0182 -0.0182 -0.0182 -0.0182 0.0182 ...
 $ s_1:List of 5
  ..$ : num [1:10] 0.634 -0.634 0.634 -0.634 0.634 ...
  ..$ : num [1:10] -0.00407 0.00407 -0.00407 0.00407 -0.00407 ...
  ..$ : num [1:10] 0.00309 -0.00309 0.00309 -0.00309 0.00309 ...
  ..$ : num [1:10] -0.00222 0.00222 -0.00222 0.00222 -0.00222 ...
  ..$ : num [1:10] 0.00744 -0.00744 0.00744 -0.00744 0.00744 ...
 $ s_2:List of 5
  ..$ : num [1:10] 0.277 -0.277 0.277 -0.277 -0.277 ...
  ..$ : num [1:10] 0.000193 -0.000193 0.000193 -0.000193 -0.000193 ...
  ..$ : num [1:10] 8.14e-05 -8.14e-05 8.14e-05 -8.14e-05 -8.14e-05 ...
  ..$ : num [1:10] 4.31e-05 -4.31e-05 4.31e-05 -4.31e-05 -4.31e-05 ...
  ..$ : num [1:10] -0.000135 0.000135 -0.000135 0.000135 0.000135 ...
 $ s_3:List of 5
  ..$ : num [1:10] 0.121 -0.121 0.121 -0.121 0.121 ...
  ..$ : num [1:10] -9.17e-06 9.17e-06 -9.17e-06 9.17e-06 -9.17e-06 ...
  ..$ : num [1:10] 2.15e-06 -2.15e-06 2.15e-06 -2.15e-06 2.15e-06 ...
  ..$ : num [1:10] -8.4e-07 8.4e-07 -8.4e-07 8.4e-07 -8.4e-07 ...
  ..$ : num [1:10] 2.46e-06 -2.46e-06 2.46e-06 -2.46e-06 2.46e-06 ...
 $ s_4:List of 5
  ..$ : num [1:10] 0.176 0.176 0.176 0.176 -0.176 ...
  ..$ : num [1:10] -7.86e-07 -7.86e-07 -7.86e-07 -7.86e-07 7.86e-07 ...
  ..$ : num [1:10] 2.51e-07 2.51e-07 2.51e-07 2.51e-07 -2.51e-07 ...
  ..$ : num [1:10] -9.56e-08 -9.56e-08 -9.56e-08 -9.56e-08 9.56e-08 ...
  ..$ : num [1:10] -1.01e-06 -1.01e-06 -1.01e-06 -1.01e-06 1.01e-06 ...
 $ s_5:List of 5
  ..$ : num [1:10] 0.2972 0.0547 0.2972 0.0547 -0.0547 ...
  ..$ : num [1:10] -9.96e-06 8.38e-06 -9.96e-06 8.38e-06 -8.38e-06 ...
  ..$ : num [1:10] 2.40e-06 -1.89e-06 2.40e-06 -1.89e-06 1.89e-06 ...
  ..$ : num [1:10] -9.35e-07 7.44e-07 -9.35e-07 7.44e-07 -7.44e-07 ...
  ..$ : num [1:10] 1.45e-06 -3.46e-06 1.45e-06 -3.46e-06 3.46e-06 ...

4 Run large simulation

Increase the parameters to realistic values and run the simulation.

# Set global parameters
v_dim <- 1e4 # vector dimension of state variables
n_step <- 1e3 # number of time steps in simulation

# Generate input variables (random)
input_df <- tibble::tibble(
  i_x = rnorm(n = n_step),
  i_y = rnorm(n = n_step)
)

# Make the initial state
init_state <- mk_init_state(v_dim, seed = 42)

# Run the simulation, returning the history as a data frame
# Time the execution with tictoc::
tictoc::tic()
history_df <- input_df %>% run_simulation(init_state, update_state)
tictoc::toc()

4.879 sec elapsed

# Take a look at the history
history_df

# A tibble: 1,000 × 11
       t    i_x     i_y    s_x    s_y s_0     s_1    s_2    s_3    s_4    s_5   
   <int>  <dbl>   <dbl>  <dbl>  <dbl> <list>  <list> <list> <list> <list> <list>
 1     1  1.30  -0.298   1.30  -0.298 <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 2     2  2.29  -0.284   3.59  -0.582 <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 3     3 -1.39   0.870   2.20   0.287 <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 4     4 -0.279 -0.544   1.92  -0.257 <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 5     5 -0.133  0.629   1.79   0.372 <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 6     6  0.636 -1.42    2.43  -1.05  <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 7     7 -0.284 -1.23    2.14  -2.28  <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 8     8 -2.66  -1.67   -0.514 -3.95  <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
 9     9 -2.44   0.0844 -2.95  -3.87  <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
10    10  1.32  -0.206  -1.63  -4.07  <dbl [… <dbl … <dbl … <dbl … <dbl … <dbl …
# … with 990 more rows

# How much RAM does the history use?
format(object.size(history_df), standard = "IEC", units = "GiB")

[1] "0.4 GiB"

That’s adequately fast.

The size is as expected.

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 21.04

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
 [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
 [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
[1] Matrix_1.3-4       DiagrammeR_1.0.6.1 tictoc_1.0.1       ggplot2_3.3.5     
[5] dplyr_1.0.7        purrr_0.3.4        tibble_3.1.3      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.7         RColorBrewer_1.1-2 pillar_1.6.2       compiler_4.1.0    
 [5] later_1.2.0        git2r_0.28.0       workflowr_1.6.2    tools_4.1.0       
 [9] digest_0.6.27      lattice_0.20-44    jsonlite_1.7.2     evaluate_0.14     
[13] lifecycle_1.0.0    gtable_0.3.0       pkgconfig_2.0.3    rlang_0.4.11      
[17] rstudioapi_0.13    cli_3.0.1          yaml_2.2.1         xfun_0.24         
[21] withr_2.4.2        stringr_1.4.0      knitr_1.33         htmlwidgets_1.5.3 
[25] generics_0.1.0     fs_1.5.0           vctrs_0.3.8        rprojroot_2.0.2   
[29] grid_4.1.0         tidyselect_1.1.1   here_1.0.1         glue_1.4.2        
[33] R6_2.5.0           fansi_0.5.0        rmarkdown_2.9      bookdown_0.22     
[37] magrittr_2.0.1     whisker_0.4        scales_1.1.1       promises_1.2.0.1  
[41] ellipsis_0.3.2     htmltools_0.5.1.1  colorspace_2.0-2   renv_0.14.0       
[45] httpuv_1.6.1       utf8_1.2.2         stringi_1.7.3      visNetwork_2.0.9  
[49] munsell_0.5.0      crayon_1.4.1

Simulation Framework

Ross Gayler

2021-07-19

1 `run_simulation()`

2 Set up example

2.1 Set global parameters

2.2 Input variables

2.3 Initial state

2.4 State update function

3 Run small simulation

4 Run large simulation

Simulation Framework

Ross Gayler

2021-07-19

1 run_simulation()

2 Set up example

2.1 Set global parameters

2.2 Input variables

2.3 Initial state

2.4 State update function

3 Run small simulation

4 Run large simulation

1 `run_simulation()`