An update to announce the migration to DifferentiationInterface.jl. Many thanks to @gdalle for all the help!
In addition to a simpler and more maintainable codebase, it allows to switch the differentiation backend for gradients and Jacobians inside NonLinMPC, MovingHorizonEstimator, linearize and ExtendedKalmanFilter. Sparse Jacobians are also supported with AutoSparse. Dense ForwardDiff.jl computation are used everywhere by default, except for the MultipleShooting transcription that uses sparse computations. Note that for small problems like the inverted pendulum with H_p=20 and H_c=2, dense Jacobians may be slightly faster than sparse matrices, even with a MultipleShooting transcription. At least, that’s what I benchmarked for this case study.
Note that the implementation rely on the Cache feature of DI.jl to reduce the allocations, and some backend does not support it for now.
The change log since my last post is:
added: migration to DifferentiationInterface.jl
added: new gradient and jacobian keyword arguments for NonLinMPC
added: new gradient and jacobian keyword arguments for MovingHorizonEstimator
added: new jacobian keyword argument for NonLinModel (for linearization)
added: new jacobian keyword argument for ExtendedKalmanFilter
added: ExtendedKalmanFilter is now allocation-free at runtime
changed: deprecate preparestate!(::SimModel,_,_), replaced by preparestate!(::SimModel)
debug: nonlinear inequality constraint with MultipleShooting now work as expected (custom + output + terminal constraints)
debug: x_noise argument in sim! now works as expected
doc: now using DocumenterInterLinks.jl to ease the maintenance
test: many new test with AutoFiniteDiff backend
test: new test to cover nonlinear inequality constraint with MultipleShooting corner cases
A quick update on the new stuff in the package since my last post.
First, the newest release introduces the ManualEstimator to turn off built-in state estimation and provide your own estimate. A first use case is to implement a linear MPC (with an approximate plant model, for the speed) with a nonlinear state estimator (with a high fidelity plant model, for accuracy). A second use case is using the exclusive observers from LowLevelParticleFilters.jl to estimate the state of the plant model and its disturbances.
Also, a significant performance boost for NonLinMPC and MovingHorizonEstimator was introduced in v1.6.0 by the more efficient value_and_gradient! and value_and_jacobian! of DI.jl. This is equivalent of using DiffResults.jl, but agnostic of the differentiation backend. I benchmarked about a 1.25x speed boost on the pendulum example of the manual.
Lastly, the nint_u option for the MovingHorizonEstimator was not working well because of a bug when the observation window is not filled (at the beginning). The bug was corrected in v1.6.2 (with new unit tests).
The next release will introduce custom move blocking, which is a way of specifying long control horizon H_c without increasing the number of decision variables in the optimization problem.
The changelog since my last post is:
added: ManualEstimator to turn off built-in state estimation and provide your own estimate \mathbf{\hat{x}}_{k}(k) or \mathbf{\hat{x}}_{k-1}(k)
added: slightly improve NonLinMPC performances with specialized conversion and weight matrices
added: significant performance boost of NonLinMPC and MovingHorizonEstimator using value_and_gradient!/jacobian! of DifferentiationInterface.jl instead of individual calls
added: setstate! now allows manual modifications of the estimation error covariance \mathbf{\hat{P}} (if computed by the estimator)
changed: M_Hp, N_Hc and L_Hp keyword arguments now default to Diagonal instead of diagm matrices for all PredictiveController constructors
changed: moved lastu0 inside PredictiveController objects
removed: DiffCaches in RungeKutta solver
debug: force update of gradient/jacobian in MovingHorzionEstimator when window not filled
debug: remove .data in KalmanFilter matrix products
debug: do not call jacobian! if nd==0 in linearize
debug: no more noisy @warn about DiffCache chunk size
test: new tests for ManualEstimator
test: added allocations tests for types that are known to be allocation-free (SKIP THEM FOR NOW)
test: adapt tests for the new automatically balancing minreal function
Here’s a quick update on the new stuff in the package.
First v1.8.0 introduces custom move blocking patterns. This parameter changes the duration in time steps of each move blocks in the profile of the manipulated inputs \mathbf{u}, which is more general than classical control horizons H_c. See the H_c argument in e.g. LinMPC for detail. Note that an intrinsic side-effect of custom move blocking is a reduced accuracy for warm-starting the optimizer.
Also, the last two versions incorporate some debugging and performance improvements for MovingHorizonEstimator. It is common that covariance matrices are strictly diagonal, especially with trial-and-error tuning. The code now preserves Diagonal types to accelerate computations in the objective function. Also, problematic termination status at optimization will not crash anymore (it will emits an @error and return the open-loop estimation).
For now, I continuously monitor the the performance of OSQP.jl, DAQP.jl, Ipopt.jl and MadNLP.jl. Note that MadNLP.jl is generally 2-3 times faster than Ipopt.jl, but it does not work well for tightly constrained problem (e.g.: multiple shooting). It’s probably related to its L-BFGS approximation. I would be curious to see the performance with exact Hessians, but it’s not supported by ModelPredictiveControl.jl for now. For differentiation, only dense and sparse ForwardDiff.jl is monitored, but I may add other backend in the future.
BTW @MilesCranmer, is there a way to preserve the order of definition in SUITE in the outputted table? Or maybe specify the order manually? Some benchmark are interlinked and I would like them to be near each other.
I will register the version soon.
Change log:
added: move blocking feature in LinMPC, ExplicitMPC and NonLinMPC (see Hc argument)
added: new KalmanCovariance parametric struct to handle covariance sparsity efficiently
added: dispatch on repeatdiag to preserve Diagonnals
added: in-place and allocation-free inv! for MHE covariance matrices
added: new benchmark suite that runs automatically on each PR with AirspeedVelocity.jl
changed: store P̃Δu, P̃u and Tu conversion matrices as SparseMatrixCSC
debug: support Hermitian weights in PredictiveController
debug: use dummy b vector in MHE setmodel! to avoid ±Inf values
debug: MovingHorizonEstimator no longer crash with error termination flag
test: verify setmodel! with He>1 for MHE
test: new integration with ManualEstimator and MPCs
test: improve coverage with error termination status
I just stumbled upon this post today and spent the whole 17 minutes reading through the developments of the package. I just wanted to say that it was a pleasure tracking all your development and looking forward to what comes up afterwards!
My only 2 cents, since you were already looking into coming up with other integrators natively beyond RK4, orthogonal collocation using Gauss-Legendre quadrature is not much more complicated and might earn you some extra stability at not a significant cost.
Another thing that might be of interest is including robust (N)MPC via scenario trees, which has been done here https://www.do-mpc.com
Thanks for the feedback, really appreciated! Yes my next step is supporting collocation as a new transcription method. I will probably implement trapezoidal integration first for its simplicity, but Gauss-Legendre will come after.
For robust MPC, it is not planned in the short term but maybe if I have time in the mid-long term.
The new release introduces a new TrapezoidalCollocation to handle moderately stiff NonLinModel!
This is presumably the simplest form of all the direct collocation methods. It internally uses the implicit trapezoidal rule, which is efficiently solved in the nonlinear equality constraints of the NLP. This is the transcription method in MATLAB’s nlmpc when the dynamics are continuous by the way. It currently assumes piecewise constant manipulated inputs \mathbf{u} (a.k.a. ZOH, like all the other integrators in the package), but I will add linear interpolation soon. The number of decision variables is identical to MultipleShooting, so the computational cost are similar. As a matter of fact, the performances are even a little bit better than MultipleShooting with RungeKutta(4):
I also did some tests with the new experimental Ipopt._VectorNonlinearOracle. This is very promising since I benchmarked a 2.5-fold speedup on NonLinMPC with MultipleShooting and custom constraints. It will also allow exact Hessians for the Lagragian function, which in turn can improve the performance if the second-order differentiation is not too heavy (likely the case with AutoSparse). I will wait for an official release from JuMP.jl, but I’m super excited to incorporate these changes on the main branch!
This release helps me a lot in understanding collocation methods. I should be able to add orthogonal collocation with Gauss-Legendre in the next releases.
The changelog since my last post:
added: TrapezoidalCollocation method for continuous nonlinear systems
added: dependabot.yml file to help with the CI dependencies
changed: new fields with the continuous state-space function model.f! and model.h! in NonLinModels (instead of model.solver_f! and model.solver_h!)
test: new test with TrapezoidalCollocation and Ipopt
bench: new case studies with TrapezoidalCollocation on the pendulum
First, the version 1.9.1 introduced first-order interpolation on the manipulated inputs \mathbf{u} with TrapezodialCollocation(1), saving one call of the state derivative function \mathbf{f} per time step (without the new multi-threading option). Note that the other ODE solvers in the package assume a zero-order hold, so there will be a slight mismatch between the controller model and a plant simulated by these integrators.
The 1.10 version added the support for a constant arrival covariance \mathbf{\bar{P}} in the MovingHorizonEstimator by constructing it with covestim=SteadyKalmanFilter(model). Doing so will fix it at the steady-state covariance \mathbf{\hat{P}}(\infty), a common choice in the litterature, but its also possible to customize the matrix.
Lastly, the version 1.11 introduces multi-threading supports for the MultipleShooting and TrapezoidalCollocation transcription, in the case of NonLinModel (linear plant model will be multithreaded if the quadratic optimizer is multithreaded). There are new f_threads and h_threads keyword arguments to activate multi-threading when calling the state \mathbf{f} and output \mathbf{h} functions in the optimization, respectively. These options are interesting if they are expensive to evaluate. This can be the case when calling a fancy implicit solver from OrdinaryDiffEq.jl in the \mathbf{f} function, for example.
I will register the release soon.
Change log:
added: multithreading with MultipleShooting and TrapezoidalCollocation
added: covestim in MovingHorizonEstimator now supports SteadyKalmanFilter for a constant arrival covariance \mathbf{\bar{P}}
added: clearer pretty-print of objects with a tree structure
added: print differentiation backends in NonLinModel, ExtendedKalmanFilter, MovingHorizonEstimator and NonLinMPC
added: piecewise linear manipulated inputs \mathbf{u} with TrapezoidalCollocation(1)
debug: setmodel! now works with LinMPC based on MultipleShooting
debug: linearize! on NonLinModel with non-zero operating points now works
test: multithreading with MultipleShooting and TrapezoidalCollocation
test: new test for setmodel! with MultipleShooting
bench: case studies with multithreading on the pendulum
doc: updated the pendulum tutorial to ModelingToolkit.jl v10 in the manual
doc: using MultipleShooting in the manual on the unstable pendulum system
doc: describing the internal of the various prediction and transcription methods
The latest update now supports the new MathOptInterface.VectorNonlinearOracle for the nonlinear constraints. There is a new oracle keyword argument in NonLinMPC and MovingHorizonEstimator to activate this feature. It is activated only for Ipopt.jl by default, since many optimizers do not support it for now. It’s about 2.5 times faster under nonlinear constraints for the pendulum case study
It’s obviously not faster with single shooting and w/o custom constraint (since there are no NL constraints), but it is worth noting that the multiple shooting is now faster than the single shooting transcription, thus the benchmarks are now in phase with the literature!
Next step is adding the support for exact Hessians.
Change log:
added: oracle=true option in NonLinMPC to significantly improve performances under nonlinear constraints (default with Ipopt.jl)
added: oracle=true option in MovingHorizonEstimator to significantly improve performances under nonlinear constraints (default with Ipopt.jl)
added: ill-conditioned QP warning at controller construction
added: progress keyword argument in sim! to disable progress logs
removed: less useful workload in precompile.jl to accelerate precompilation
doc: minor precision in strictly proper system assumption
It’s been a long time since my last post, so let’s post a little update on the new features and bugfixes.
The support for exact Hessian was added in v1.13, and the following patches introduced performance tweaks on sparse matrix colouring. Exact Hessians can be activated in NonLinMPC and MovingHorizonEstimator at construction with the option hessian=true. Otherwise, the internal quasi-Newton method of the optimizer will be used, as it was the case before. It is also now possible to inspect the gradients, Jacobians and Hessians at the optimum using the getinfo! method.
The last two patches included one bugfix on the default softness parameters of the terminal constraints \mathbf{c_{\hat{x}_{min}}} and \mathbf{c_{\hat{x}_{max}}}. They can lead to infeasible problems so they should be soft by default, and they were hard by default because of a bug that is now fixed.
The last bugfix was related to custom move blocking patterns using LinModel, affecting both LinMPC and NonLinMPC controllers. The solution of the optimization problem was simply dead wrong, because of a mistake at the construction of the prediction equations and matrices. Turns out that the math are way more complex than I first thought. I needed to revise my notation to rigorously keep track of all the various indices in the block matrices. The new notation is documented in init_predmat and init_defectmat, in the extended help sections. BTW, it took me around 50 hours to derive the math correctly, implement it, and add new tests with custom move blocking!
I hit the last bug while I was working on an interface with LinearMPC.jl, to enable lightweight C-code generation for embedded device. The next release should introduce this new interface!
I will register the latest release soon. Happy new year to everyone !
The change log since my last post:
added: hessian argument in NonLinMPC for exact Hessian of Lagrangian
added: hessian argument in MovingHorizonEstimator for exact Hessian of Lagrangian
added: trunc_cov function for MovingHorizonEstimator (to optimize perf. of dot in objective function)
added: lastu keyword argument to moveinput!
added: derivatives in getinfo dictionnary for NonLinMPC and MovingHorizonEstimator
added: optimal nonlinear constraint vectors in getinfo
changed: new default sparse backend for jacobian and hessian in NonLinMPC and MovingHorizonEstimator
debug: prepare AD with rand constant contexts
debug: only fills non-zeros in ∇²f of @operator
debug: default terminal constraint softness parameters c_x̂min/max equal to 1
debug: correctly handle custom move blocking for LinModel
The new release introduces an interface with LinearMPC.jl for the LinMPC controllers, enabling lightweight C code generation for embedded platforms. These features are implemented as a package extension that will be loaded if LinearMPC is present in the current active environment. Many thanks to @darnstrom for all the help during the development!
To do so, one first creates a ModelPredictiveControl.LinMPC object, convert it to LinearMPC.MPC object, and call codegen on the latter:
Note that not all features of ModelPredictiveControl.jl are supported, see the documentation for details on this aspect. But the reverse is also true: LinMPC.jl offers exclusive features such as binary manipulated inputs and constrained explicit MPC.
I will probably work on new collocation methods for the next releases.
Change log:
added: C code generation via LinearMPC.jl
doc: various improvement in the manual and function sections
test: new closed-loop tests with LinMPC converted to LinearMPC.MPC
Very cool! Have you done any profiling of flash/ram requirements for the generated C code?
For a 500 variable model I currently need about 300 kB flash / RAM on an STM32 where the C code is generated by cvxpygen (and the MPC is defined in cvxpy), on float64.
No I did not, maybe @darnstrom has some information on this aspect ? I do have an old STM32 lying around somewhere at home, that could be an interesting project for me! It’s been a long time since I’ve worked on this platform. What do you use for the IDE/compiler/uploader @langestefan (preferably on linux) ?
Do those 500 decision variables include states? LinearMPC.jl eliminates the states, so only the controls are decision variables, but from my experience this is not so easy to do with the current version of CVXPY/CVXPYgen (which can lead to an inflated memory footprint.)
From my experience, the memory footprint and execution time from the code generated by LinearMPC.jl can be significantly faster than code from CVXPYgen for MPC problems, but producing concrete results to strengthen this is still a TODO (I just have some MPC applications I have run myself where I have observed this) I have plans to do a more formal comparison of the generated code with, for example, CVXPYgen, acados, and other packages that support code generation)
But there is a trade-off here right? Eliminating the state variables changes the structure of your problem matrices. I found that with a good sparse QP solver extra variables are not necessarily a problem. The nice thing about cvxpy+cvxpygen is that you can try a bunch of solvers and the codegen works all the same.
For the particular application I am working on robustness of the solver and code size are the most important criteria. The dynamics are slow, so solve time is not important.
I am very curious now to see how LinearMPC.jl + ModelPredictiveControl.jl compares to my current setup
Maybe not a problem for the solution time, but they still lead to large (even if sparse) matrices that must be stored? Take the extreme case of a 1 horizon problem with 1 control input and 500-dimensional state. The QP in LinearMPC would be 1x1, while a sparse problem with state among the decision variables will have at least 500 non-zeros in the hessian?
In terms of numerical robustness, the original paper of DAQP showed that it performed better when the condition number of the Hessian is fairly large (> 10^{10}) and when the prediction horizon is not too high (as it should be in most cases). It also aligns with my experience, although sparse solvers like OSQP can be “robustified” by introducing the state in the decision variables (a.k.a. multiple shoothing, I think all case studies were in a single shooting form in the paper). See here in which I compare:
OSQP.jl in single shooting
DAQP.jl in single shooting
OSQP.jl in multiple shooting
The results of 2. and 3. are very similar on this simple case study. I expect that the robustness will be similar to your current setup.
Not sure if I am missing something here, but does ModelPredictiveControl.jl (and LinearMPC.jl) not support time-varying constraints on inputs and states?
