diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index cf4d5bc..99fa545 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-09-09T16:05:58","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-09-09T16:17:25","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/examples/PG_OCP_generic_basis_functions/index.html b/dev/examples/PG_OCP_generic_basis_functions/index.html index fea525c..a311470 100644 --- a/dev/examples/PG_OCP_generic_basis_functions/index.html +++ b/dev/examples/PG_OCP_generic_basis_functions/index.html @@ -274,4 +274,4 @@ end # Plot predictions. -plot_predictions(y_opt, y_sys; plot_percentiles=false, y_min=y_min, y_max=y_max) +plot_predictions(y_opt, y_sys; plot_percentiles=false, y_min=y_min, y_max=y_max) diff --git a/dev/examples/PG_OCP_known_basis_functions/index.html b/dev/examples/PG_OCP_known_basis_functions/index.html index 47241d4..c03e411 100644 --- a/dev/examples/PG_OCP_known_basis_functions/index.html +++ b/dev/examples/PG_OCP_known_basis_functions/index.html @@ -204,4 +204,4 @@ end # Plot predictions. -plot_predictions(y_opt, y_sys; plot_percentiles=false, y_min=y_min, y_max=y_max) +plot_predictions(y_opt, y_sys; plot_percentiles=false, y_min=y_min, y_max=y_max) diff --git a/dev/examples/autocorrelation/index.html b/dev/examples/autocorrelation/index.html index 016b9c8..fb5c129 100644 --- a/dev/examples/autocorrelation/index.html +++ b/dev/examples/autocorrelation/index.html @@ -63,4 +63,4 @@ y_test = y[:, T+1:end]

Infer model

Run the particle Gibbs sampler to jointly estimate the model parameters and the latent state trajectory.

PG_samples = particle_Gibbs(u_training, y_training, K, K_b, k_d, N, phi, Lambda_Q, ell_Q, Q_init, V, A_init, x_init_dist, g, R)
 
-time_sampling = time() - sampling_timer

Plot autocorrelation

Finally, plot the autocorrelation of the obtained samples.

plot_autocorrelation(PG_samples; max_lag=100)
+time_sampling = time() - sampling_timer

Plot autocorrelation

Finally, plot the autocorrelation of the obtained samples.

plot_autocorrelation(PG_samples; max_lag=100)
diff --git a/dev/index.html b/dev/index.html index 8e10cad..a2fcf96 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Introduction ·
GitHub

Introduction

PGopt is a software for determining optimal input trajectories with probabilistic performance and constraint satisfaction guarantees for unknown systems with latent states based on input-output measurements. In order to quantify uncertainties, which is crucial for deriving formal guarantees, a Bayesian approach is employed and a prior over the unknown dynamics and the system trajectory is formulated in state-space representation. Since for practical applicability, the prior must be updated based on input-output measurements, but the corresponding posterior distribution is analytically intractable, particle Gibbs (PG) sampling is utilized to draw samples from this distribution. Based on these samples, a scenario optimal control problem (OCP) is formulated and probabilistic performance and constraint satisfaction guarantees are inferred via a greedy constraint removal.

The approach is explained in the paper "Learning-Based Optimal Control with Performance Guarantees for Unknown Systems with Latent States", available as a preprint on arXiv.

Versions

This document describes the Julia implementation of PGopt. In this version, the solvers Altro and IPOPT can be employed to solve the optimal control problem (OCP).

The results presented in the paper were generated with the solver Altro. Please note that this implementation has several limitations: only cost functions of the form $J_H=\sum\nolimits_{t=0}^H \frac{1}{2} u_t R u_t$, measurement functions of the form $y=x_{1:n_y}$, and output constraints of the form $y_\mathrm{min} \leq y \leq y_\mathrm{max}$ are supported.

The solver IPOPT is more general than Altro, and this implementation allows arbitrary cost functions $J_H(u_{1:H},x_{1:H},y_{1:H})$, measurement functions $y=g(x,u)$, and constraints $h(u_{1:H},x_{1:H},y_{1:H})$. However, using IPOPT together with the proprietary HSL Linear Solvers for Julia (HSL_jll.jl) is recommended. A license (free to academics) is required.

Besides the Julia implementation, there is also a MATLAB implementation that utilizes CasADi and IPOPT. Further information can be found here.

Installation

PGopt can be installed using the Julia package manager. Start a Pkg REPL (press ] in a Julia REPL), and install PGopt via

pkg> add https://github.com/TUM-ITR/PGopt:Julia

Alternatively, to inspect the source code more easily, download the source code from GitHub. Navigate to the folder PGopt/Julia, start a Pkg REPL (press ] in a Julia REPL), and install the dependencies via

pkg>activate . 
-pkg>instantiate

You can then execute the examples in the folder PGopt/Julia/examples.

+pkg>instantiate

You can then execute the examples in the folder PGopt/Julia/examples.

diff --git a/dev/list_of_functions/index.html b/dev/list_of_functions/index.html index 834a643..7db8f53 100644 --- a/dev/list_of_functions/index.html +++ b/dev/list_of_functions/index.html @@ -1,30 +1,30 @@ List of fuctions ·
GitHub

List of fuctions

PGopt.particle_Gibbs — Method
particle_Gibbs(u_training, y_training, K, K_b, k_d, N, phi::Function, Lambda_Q, ell_Q, Q_init, V, A_init, x_init_dist, g, R; x_prim=nothing)

Run particle Gibbs sampler with ancestor sampling to obtain samples $\{A, Q, x_{T:-1}\}^{[1:K]}$ from the joint parameter and state posterior distribution $p(A, Q, x_{T:-1} \mid \mathbb{D}=\{u_{T:-1}, y_{T:-1}\})$.

Arguments

  • u_training: training input trajectory
  • y_training: training output trajectory
  • K: number of models/scenarios to be sampled
  • K_b: length of the burn in period
  • k_d: number of models/scenarios to be skipped to decrease correlation (thinning)
  • N: number of particles
  • phi: basis functions
  • Lambda_Q: scale matrix of IW prior on $Q$
  • ell_Q: degrees of freedom of IW prior on $Q$
  • Q_init: initial value of $Q$
  • V: left covariance matrix of MN prior on $A$
  • A_init: initial value of $A$
  • x_init_dist: distribution of the initial state
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise
  • x_prim: prespecified trajectory for the first iteration

This function is based on the papers

A. Svensson and T. B. Schön, “A flexible state–space model for learning nonlinear dynamical systems,” Automatica, vol. 80, pp. 189–199, 2017.
 
-F. Lindsten, T. B. Schön, and M. Jordan, “Ancestor sampling for particle Gibbs,” Advances in Neural Information Processing Systems, vol. 25, 2012.

and the code provided in the supplementary material.

source
PGopt.solve_PG_OCP_Altro — Method
solve_PG_OCP_Altro(PG_samples::Vector{PG_sample}, phi::Function, R, H, u_min, u_max, y_min, y_max, R_cost_diag; x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, active_constraints=nothing, opts=nothing, print_progress=true)

Solve the optimal control problem of the following form using Altro.jl:

$\min \sum_{t=0}^{H} \frac{1}{2} u_t \operatorname{diag}(R_{\mathrm{cost}}) u_t$

subject to:

\[\begin{aligned} +F. Lindsten, T. B. Schön, and M. Jordan, “Ancestor sampling for particle Gibbs,” Advances in Neural Information Processing Systems, vol. 25, 2012.

and the code provided in the supplementary material.

source
PGopt.solve_PG_OCP_Altro — Method
solve_PG_OCP_Altro(PG_samples::Vector{PG_sample}, phi::Function, R, H, u_min, u_max, y_min, y_max, R_cost_diag; x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, active_constraints=nothing, opts=nothing, print_progress=true)

Solve the optimal control problem of the following form using Altro.jl:

$\min \sum_{t=0}^{H} \frac{1}{2} u_t \operatorname{diag}(R_{\mathrm{cost}}) u_t$

subject to:

\[\begin{aligned} \forall k, \forall t \\ x_{t+1}^{[k]} &= f_{\theta^{[k]}}(x_t^{[k]}, u_t) + v_t^{[k]}, \\ x_{t, 1:n_y}^{[k]} &\geq y_{\mathrm{min},\ t} - e_t^{[k]}, \\ x_{t, 1:n_y}^{[k]} &\leq y_{\mathrm{max},\ t} - e_t^{[k]}, \\ u_t &\geq u_{\mathrm{min},\ t}, \\ u_t &\leq u_{\mathrm{max},\ t}. -\end{aligned}\]

Note that the output constraints imply the measurement function $y_t^{[k]} = x_{t, 1:n_y}^{[k]}$. Further note that the states of the individual models ($x^{[1:K]}$) are combined in the vector x_vec of dimension K * n_x.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • u_min: array of dimension 1 x 1 or 1 x H containing the minimum control input for all timesteps
  • u_max: array of dimension 1 x 1 or 1 x H containing the maximum control input for all timesteps
  • y_min: array of dimension 1 x 1 or 1 x H containing the minimum system output for all timesteps
  • y_max: array of dimension 1 x 1 or 1 x H containing the maximum system output for all timesteps
  • R_cost_diag: parameter of the diagonal quadratic cost function
  • x_vec_0: vector with $K n_x$ elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • active_constraints: vector containing the indices of the models, for which the output constraints are active - if not provided, the output constraints are considered for all models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Altro_greedy_guarantees — Method
solve_PG_OCP_Altro_greedy_guarantees(PG_samples::Vector{PG_sample}, phi::Function, R, H, u_min, u_max, y_min, y_max, R_cost_diag, β; x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, opts=nothing, print_progress=true)

Solve the following optimal control problem using Altro.jl and determine a support sub-sample with cardinality $s$ via a greedy constraint removal. Based on the cardinality $s$, a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory $u_{0:H}$ is applied to the unknown system is calculated (i.e., $1-\epsilon$ is determined).

$\min \sum_{t=0}^{H} \frac{1}{2} u_t \operatorname{diag}(R_{\mathrm{cost}}) u_t$

subject to:

\[\begin{aligned} +\end{aligned}\]

Note that the output constraints imply the measurement function $y_t^{[k]} = x_{t, 1:n_y}^{[k]}$. Further note that the states of the individual models ($x^{[1:K]}$) are combined in the vector x_vec of dimension K * n_x.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • u_min: array of dimension 1 x 1 or 1 x H containing the minimum control input for all timesteps
  • u_max: array of dimension 1 x 1 or 1 x H containing the maximum control input for all timesteps
  • y_min: array of dimension 1 x 1 or 1 x H containing the minimum system output for all timesteps
  • y_max: array of dimension 1 x 1 or 1 x H containing the maximum system output for all timesteps
  • R_cost_diag: parameter of the diagonal quadratic cost function
  • x_vec_0: vector with $K n_x$ elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • active_constraints: vector containing the indices of the models, for which the output constraints are active - if not provided, the output constraints are considered for all models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Altro_greedy_guarantees — Method
solve_PG_OCP_Altro_greedy_guarantees(PG_samples::Vector{PG_sample}, phi::Function, R, H, u_min, u_max, y_min, y_max, R_cost_diag, β; x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, opts=nothing, print_progress=true)

Solve the following optimal control problem using Altro.jl and determine a support sub-sample with cardinality $s$ via a greedy constraint removal. Based on the cardinality $s$, a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory $u_{0:H}$ is applied to the unknown system is calculated (i.e., $1-\epsilon$ is determined).

$\min \sum_{t=0}^{H} \frac{1}{2} u_t \operatorname{diag}(R_{\mathrm{cost}}) u_t$

subject to:

\[\begin{aligned} \forall k, \forall t \\ x_{t+1}^{[k]} &= f_{\theta^{[k]}}(x_t^{[k]}, u_t) + v_t^{[k]}, \\ x_{t, 1:n_y}^{[k]} &\geq y_{\mathrm{min},\ t} - e_t^{[k]}, \\ x_{t, 1:n_y}^{[k]} &\leq y_{\mathrm{max},\ t} - e_t^{[k]}, \\ u_t &\geq u_{\mathrm{min},\ t}, \\ u_t &\leq u_{\mathrm{max},\ t}. -\end{aligned}\]

Note that the output constraints imply the measurement function $y_t^{[k]} = x_{t, 1:n_y}^{[k]}$. Further note that the states of the individual models ($x^{[1:K]}$) are combined in the vector x_vec of dimension K * n_x.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • u_min: array of dimension 1 x 1 or 1 x H containing the minimum control input for all timesteps
  • u_max: array of dimension 1 x 1 or 1 x H containing the maximum control input for all timesteps
  • y_min: array of dimension 1 x 1 or 1 x H containing the minimum system output for all timesteps
  • y_max: array of dimension 1 x 1 or 1 x H containing the maximum system output for all timesteps
  • R_cost_diag: parameter of the diagonal quadratic cost function
  • β: confidence parameter
  • x_vec_0: vector with K*n_x elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Ipopt — Method
solve_PG_OCP_Ipopt(PG_samples::Vector{PG_sample}, phi::Function, g::Function, R, H, J, h_scenario, h_u; J_u=false, x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, solver_opts=nothing, print_progress=true)

Solve the optimal control problem of the following form using Ipopt:

$\min_{u_{0:H},\; \overline{J_H}} \overline{J_H}$

subject to:

\[\begin{aligned} +\end{aligned}\]

Note that the output constraints imply the measurement function $y_t^{[k]} = x_{t, 1:n_y}^{[k]}$. Further note that the states of the individual models ($x^{[1:K]}$) are combined in the vector x_vec of dimension K * n_x.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • u_min: array of dimension 1 x 1 or 1 x H containing the minimum control input for all timesteps
  • u_max: array of dimension 1 x 1 or 1 x H containing the maximum control input for all timesteps
  • y_min: array of dimension 1 x 1 or 1 x H containing the minimum system output for all timesteps
  • y_max: array of dimension 1 x 1 or 1 x H containing the maximum system output for all timesteps
  • R_cost_diag: parameter of the diagonal quadratic cost function
  • β: confidence parameter
  • x_vec_0: vector with K*n_x elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Ipopt — Method
solve_PG_OCP_Ipopt(PG_samples::Vector{PG_sample}, phi::Function, g::Function, R, H, J, h_scenario, h_u; J_u=false, x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, solver_opts=nothing, print_progress=true)

Solve the optimal control problem of the following form using Ipopt:

$\min_{u_{0:H},\; \overline{J_H}} \overline{J_H}$

subject to:

\[\begin{aligned} \forall k, &\forall t \\ x_{t+1}^{[k]} &= f_{\theta^{[k]}}(x_t^{[k]}, u_t) + v_t^{[k]}, \\ y_{t}^{[k]} &= g_{\theta^{[k]}}(x_t^{[k]}, u_t) + w_t^{[k]}, \\ J_H^{[k]} &= J_H(u_{0:H}, x_{0:H}^{[k]}, y_{0:H}^{[k]}) \leq \overline{J_H}, \\ h(&u_{0:H},x_{0:H}^{[k]},y_{0:H}^{[k]}) \leq 0. -\end{aligned}\]

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • J: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) (or $u_{1:H}$ if J_u is set true) that returns the cost to be minimized
  • h_scenario: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) that returns the constraint vector belonging to a scenario; a feasible solution must satisfy $h_{\mathrm{scenario}} \leq 0$ for all scenarios.
  • h_u: function with input argument $u_{1:H}$ that returns the constraint vector for the control inputs; a feasible solution satisfy $h_u \leq 0$.
  • J_u: set to true if cost depends only on inputs `u_{1:H} - this accelerates the optimization
  • x_vec_0: vector with K * n_x elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Ipopt_greedy_guarantees — Method
solve_PG_OCP_Ipopt_greedy_guarantees(PG_samples::Vector{PG_sample}, phi::Function, g::Function, R, H, J, h_scenario, h_u, β; J_u=false, x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, solver_opts=nothing, print_progress=true)

Solve the sample-based optimal control problem using Ipopt and determine a support sub-sample with cardinality s via a greedy constraint removal. Based on the cardinality s, a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory u_{0:H} is applied to the unknown system is calculated.

$\min_{u_{0:H},\; \overline{J_H}} \overline{J_H}$

subject to:

\[\begin{aligned} +\end{aligned}\]

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • J: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) (or $u_{1:H}$ if J_u is set true) that returns the cost to be minimized
  • h_scenario: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) that returns the constraint vector belonging to a scenario; a feasible solution must satisfy $h_{\mathrm{scenario}} \leq 0$ for all scenarios.
  • h_u: function with input argument $u_{1:H}$ that returns the constraint vector for the control inputs; a feasible solution satisfy $h_u \leq 0$.
  • J_u: set to true if cost depends only on inputs `u_{1:H} - this accelerates the optimization
  • x_vec_0: vector with K * n_x elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.solve_PG_OCP_Ipopt_greedy_guarantees — Method
solve_PG_OCP_Ipopt_greedy_guarantees(PG_samples::Vector{PG_sample}, phi::Function, g::Function, R, H, J, h_scenario, h_u, β; J_u=false, x_vec_0=nothing, v_vec=nothing, e_vec=nothing, u_init=nothing, K_pre_solve=0, solver_opts=nothing, print_progress=true)

Solve the sample-based optimal control problem using Ipopt and determine a support sub-sample with cardinality s via a greedy constraint removal. Based on the cardinality s, a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory u_{0:H} is applied to the unknown system is calculated.

$\min_{u_{0:H},\; \overline{J_H}} \overline{J_H}$

subject to:

\[\begin{aligned} \forall k, &\forall t \\ x_{t+1}^{[k]} &= f_{\theta^{[k]}}(x_t^{[k]}, u_t) + v_t^{[k]}, \\ y_{t}^{[k]} &= g_{\theta^{[k]}}(x_t^{[k]}, u_t) + w_t^{[k]}, \\ J_H^{[k]} &= J_H(u_{0:H}, x_{0:H}^{[k]}, y_{0:H}^{[k]}) \leq \overline{J_H}, \\ h(&u_{0:H},x_{0:H}^{[k]},y_{0:H}^{[k]}) \leq 0. -\end{aligned}\]

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise - only used if e_vec is not passed
  • H: horizon of the OCP
  • J: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) (or $u_{1:H}$ if J_u is set true) that returns the cost to be minimized
  • h_scenario: function with input arguments ($u_{1:H}$, $x_{1:H}$, $y_{1:H}$) that returns the constraint vector belonging to a scenario; a feasible solution must satisfy $h_{\mathrm{scenario}} \leq 0$ for all scenarios.
  • h_u: function with input argument $u_{1:H}$ that returns the constraint vector for the control inputs; a feasible solution satisfy $h_u \leq 0$.
  • β: confidence parameter
  • J_u: set to true if cost depends only on inputs `u_{1:H} - this accelerates the optimization
  • x_vec_0: vector with K * n_x elements containing the initial state of all models - if not provided, the initial states are sampled based on the PGS samples
  • v_vec: array of dimension n_x x H x K that contains the process noise for all models and all timesteps - if not provided, the noise is sampled based on the PGS samples
  • e_vec: array of dimension n_y x H x K that contains the measurement noise for all models and all timesteps - if not provided, the noise is sampled based on the provided R
  • u_init: initial guess for the optimal trajectory
  • K_pre_solve: if K_pre_solve > 0, an initial guess for the optimal trajectory is obtained by solving the OCP with only K_pre_solve < K models
  • opts: SolverOptions struct containing options of the solver
  • print_progress: if set to true, the progress is printed
source
PGopt.test_prediction — Method
test_prediction(PG_samples::Vector{PG_sample}, phi::Function, g, R, k_n, u_test, y_test)

Simulate the PGS samples forward in time and compare the predictions to the test data.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise
  • k_n: each model is simulated $k_n$ times
  • u_test: test input
  • y_test: test output
source
PGopt.plot_predictions — Method
plot_predictions(y_pred, y_test; plot_percentiles=false, y_min=nothing, y_max=nothing)

Plot the predictions and the test data.

Arguments

  • y_pred: matrix containing the output predictions
  • y_test: test output trajectory
  • plot_percentiles: if set to true, percentiles are plotted
  • y_min: min output to be plotted as constraint
  • y_max: max output to be plotted as constraint
source
PGopt.plot_autocorrelation — Method
plot_autocorrelation(PG_samples::Vector{PG_sample}; max_lag=0)

Plot the autocorrelation function (ACF) of the PG samples. This might be helpful when adjusting the thinning parameter $k_d$.

Arguments

  • PG_samples: PG samples
  • max_lag: maximum lag at which to calculate the ACF
source
PGopt.epsilon — Method
epsilon(s::Int64, K::Int64, β::Float64)

Determine the parameter $\epsilon$. $1-\epsilon$ corresponds to a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory $u_{0:H}$ is applied to the unknown system. $\epsilon$ is the unique solution over the interval $(0,1)$ of the polynomial equation in the $v$ variable:

$\binom{K}{s}(1-v)^{K-s}-\frac{\beta}{K}\sum_{m=s}^{K-1}\binom{m}{s}(1-v)^{m-s}=0$.

Arguments

  • s: cardinality of the support sub-sample
  • K: number of scenarios
  • β: confidence parameter

This function is based on the paper

S. Garatti and M. C. Campi, “Risk and complexity in scenario optimization,” Mathematical Programming, vol. 191, no. 1, pp. 243–279, 2022.

and the code provided in the appendix.

source
PGopt.MNIW_sample — Method
MNIW_sample(Phi, Psi, Sigma, V, Lambda_Q, ell_Q, T)

Sample new model parameters $\{A, Q\}$ from the conditional distribution $p(A, Q | x_{T:-1})$, which is a matrix normal inverse Wishart (MNIW) distribution.

Arguments

  • Phi: statistic; see paper below for definition
  • Psi: statistic; see paper below for definition
  • Sigma: statistic; see paper below for definition
  • V: left covariance matrix of MN prior on $A$
  • Lambda_Q: scale matrix of IW prior on $Q$
  • ell_Q: degrees of freedom of IW prior on $Q$
  • T: length of the training trajectory

This function is based on the paper

A. Svensson and T. B. Schön, “A flexible state–space model for learning nonlinear dynamical systems,” Automatica, vol. 80, pp. 189–199, 2017.

and the code provided in the supplementary material.

source
PGopt.systematic_resampling — Method
systematic_resampling(W, N)

Sample N indices according to the weights W. This function could be replaced by the function wsample from StatsBase.jl but is kept to minimize the differences between the Julia and the MATLAB implementations.

Arguments

  • W: vector containing the weights of the particles
  • N: number of indices to be sampled

This function is based on the paper

A. Svensson and T. B. Schön, “A flexible state–space model for learning nonlinear dynamical systems,” Automatica, vol. 80, pp. 189–199, 2017.

and the code provided in the supplementary material.

source
+\end{aligned}\]

Arguments

source
PGopt.test_prediction — Method
test_prediction(PG_samples::Vector{PG_sample}, phi::Function, g, R, k_n, u_test, y_test)

Simulate the PGS samples forward in time and compare the predictions to the test data.

Arguments

  • PG_samples: PG samples
  • phi: basis functions
  • g: observation function
  • R: variance of zero-mean Gaussian measurement noise
  • k_n: each model is simulated $k_n$ times
  • u_test: test input
  • y_test: test output
source
PGopt.plot_predictions — Method
plot_predictions(y_pred, y_test; plot_percentiles=false, y_min=nothing, y_max=nothing)

Plot the predictions and the test data.

Arguments

  • y_pred: matrix containing the output predictions
  • y_test: test output trajectory
  • plot_percentiles: if set to true, percentiles are plotted
  • y_min: min output to be plotted as constraint
  • y_max: max output to be plotted as constraint
source
PGopt.plot_autocorrelation — Method
plot_autocorrelation(PG_samples::Vector{PG_sample}; max_lag=0)

Plot the autocorrelation function (ACF) of the PG samples. This might be helpful when adjusting the thinning parameter $k_d$.

Arguments

  • PG_samples: PG samples
  • max_lag: maximum lag at which to calculate the ACF
source
PGopt.epsilon — Method
epsilon(s::Int64, K::Int64, β::Float64)

Determine the parameter $\epsilon$. $1-\epsilon$ corresponds to a bound on the probability that the incurred cost exceeds the worst-case cost or that the constraints are violated when the input trajectory $u_{0:H}$ is applied to the unknown system. $\epsilon$ is the unique solution over the interval $(0,1)$ of the polynomial equation in the $v$ variable:

$\binom{K}{s}(1-v)^{K-s}-\frac{\beta}{K}\sum_{m=s}^{K-1}\binom{m}{s}(1-v)^{m-s}=0$.

Arguments

  • s: cardinality of the support sub-sample
  • K: number of scenarios
  • β: confidence parameter

This function is based on the paper

S. Garatti and M. C. Campi, “Risk and complexity in scenario optimization,” Mathematical Programming, vol. 191, no. 1, pp. 243–279, 2022.

and the code provided in the appendix.

source
PGopt.MNIW_sample — Method
MNIW_sample(Phi, Psi, Sigma, V, Lambda_Q, ell_Q, T)

Sample new model parameters $\{A, Q\}$ from the conditional distribution $p(A, Q | x_{T:-1})$, which is a matrix normal inverse Wishart (MNIW) distribution.

Arguments

  • Phi: statistic; see paper below for definition
  • Psi: statistic; see paper below for definition
  • Sigma: statistic; see paper below for definition
  • V: left covariance matrix of MN prior on $A$
  • Lambda_Q: scale matrix of IW prior on $Q$
  • ell_Q: degrees of freedom of IW prior on $Q$
  • T: length of the training trajectory

This function is based on the paper

A. Svensson and T. B. Schön, “A flexible state–space model for learning nonlinear dynamical systems,” Automatica, vol. 80, pp. 189–199, 2017.

and the code provided in the supplementary material.

source
PGopt.systematic_resampling — Method
systematic_resampling(W, N)

Sample N indices according to the weights W. This function could be replaced by the function wsample from StatsBase.jl but is kept to minimize the differences between the Julia and the MATLAB implementations.

Arguments

  • W: vector containing the weights of the particles
  • N: number of indices to be sampled

This function is based on the paper

A. Svensson and T. B. Schön, “A flexible state–space model for learning nonlinear dynamical systems,” Automatica, vol. 80, pp. 189–199, 2017.

and the code provided in the supplementary material.

source
diff --git a/dev/reference/index.html b/dev/reference/index.html index 2ef4203..711560d 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -6,4 +6,4 @@ pages={90--97}, year={2024}, organization={IEEE} -} +}