This folder contains 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
The solver IPOPT is more general than Altro, and this implementation allows arbitrary cost functions HSL_jll.jl) is recommended. A license (free to academics) is required.
To execute the code, install Julia. Navigate to this folder, start a Pkg REPL (press ] in a Julia REPL), and install the dependencies via
pkg>activate .
pkg>instantiate
Then execute the example scripts in the folder examples.
Tested with Windows 11 and Julia 1.10.
The following examples can be found in the folder examples.
This script reproduces the normalized auto-correlation function plot (Figure 1) in Section V-B of the paper.
Assuming knowledge of the basis functions, samples are drawn from the posterior distribution over model parameters and latent state trajectories using the function particle_Gibbs() without thinning. Afterward, the autocorrelation is plotted using the function plot_autocorrelation().
The runtime of the script is about 10 minutes on a standard laptop.
This script reproduces the results of the optimal control approach with known basis functions (Figure 2) given in Section V-B of the paper using the solver Altro.
### Support sub-sample found
Cardinality of the support sub-sample (s): 7
Max. constraint violation probability (1-epsilon): 10.22 %
First, the algorithm and simulation parameters are defined, and training data is generated. Then, by calling the function particle_Gibbs(), samples are drawn from the posterior distribution using particle Gibbs sampling. These samples are then passed to the function solve_PG_OCP_greedy_guarantees_Altro(), which solves the scenario OCP using Altro and inferres probabilistic constraint satisfaction guarantees by greedily removing constraints and solving the corresponding reduced OCP.
The runtime of the script is about 2 hours on a standard laptop.
For the results in Table II of the paper, this script is repeated with seeds 1:100.
This script is similar to PG_OCP_known_basis_functions_Altro.jl, but the solver IPOPT is used instead of Altro. Due to the different solver, the results differ slightly from the ones presented in the paper.
This script reproduces the results of the optimal control approach with generic basis functions (Figure 3) given in Section V-C of the paper. Due to small adjustments to the code (see Erratum in the general README) and resulting numerical differences, the results from v0.2.2 and later differ slightly from those presented in the paper. If you wish to reproduce the exact results from the paper, you can pull an earlier release of the code before version v0.2.2.
In this example, the method presented in the paper "A flexible state–space model for learning nonlinear dynamical systems" is utilized to systematically derive basis functions and priors for the parameters based on a reduced-rank GP approximation. Afterward, by calling the function particle_Gibbs(), samples are drawn from the posterior distribution using particle Gibbs sampling. These samples are then passed to the function solve_PG_OCP_Altro(), which solves the scenario OCP using Altro.
The runtime of the script is about 1 hour on a standard laptop.
For the results in Table IV of the paper, this script is repeated with seeds 1:100.
This script is similar to PG_OCP_generic_basis_functions_Altro.jl, but the solver IPOPT is used instead of Altro.
The following files in the folder src contain the PGopt source code.
Defines the module PGopt, which contains the functions defined in particle_Gibbs.jl, optimal_control_Altro.jl, and optimal_control_Ipopt.jl.
This file contains the particle Gibbs sampler (function particle_Gibbs()) and several auxiliary functions that can, for example, be used to test the predictions.
| Function | Purpose |
|---|---|
| particle_Gibbs() | Run particle Gibbs sampler with ancestor sampling to obtain samples from the joint parameter and state posterior distribution. |
| test_prediction() | Simulate the PGS samples forward in time and compare the predictions to the test data. |
| plot_predictions() | Plot the predictions and the test data. |
| plot_autocorrelation() | Plot the autocorrelation function (ACF) of the PG samples. |
| systematic_resampling() | Helper function for particle_Gibbs(). Sample indices according to the weights. |
| MNIW_sample() | Helper function for particle_Gibbs(). Sample new model parameters from the distribution conditional on the previously sampled state trajectory. |
This file contains all the functions required to solve the scenario OCP and to derive the theoretical guarantees. The solver Altro is used for the optimization.
| Function | Purpose |
|---|---|
| solve_PG_OCP_Altro() | Solve the scenario optimal control problem using Altro. |
| solve_PG_OCP_Altro_greedy_guarantees() | Solve the scenario optimal control problem using Altro and determine a support sub-sample with cardinality s via a greedy constraint removal. Based on the cardinality s, calculate 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() | Determine epsilon based on the cardinality s. |
| RobotDynamics.discrete_dynamics() | Helper function for solve_PG_OCP(). |
This file contains all the functions required to solve the scenario OCP and to derive the theoretical guarantees. The solver IPOPT is used for the optimization.
| Function | Purpose |
|---|---|
| solve_PG_OCP_Ipopt() | Solve the scenario optimal control problem using IPOPT. |
| solve_PG_OCP_Ipopt_greedy_guarantees() | Solve the scenario optimal control problem using IPOPT and determine a support sub-sample with cardinality s via a greedy constraint removal. Based on the cardinality s, calculate 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. |
Detailed explanations of the individual functions can be found as comments in the file.