diff --git a/Model_Generation.m b/Model_Generation.m index ce91d34..4193cb9 100644 --- a/Model_Generation.m +++ b/Model_Generation.m @@ -21,6 +21,15 @@ run(settings.model); %% +% The code reported in this script refers to the papers reported below +% - [1] Y. Chen, M. Bruschetta, E. Picotti and A. Beghi, "MATMPC - A MATLAB +% Based Toolbox for Real-time Nonlinear Model Predictive Control," 2019 18th +% European Control Conference (ECC), 2019, pp. 3365-3370, doi: 10.23919/ECC.2019.8795788. +% - [2] Y. Chen, M. Bruschetta, D. Cuccato and A. Beghi, "An Adaptive Partial +% Sensitivity Updating Scheme for Fast Nonlinear Model Predictive Control," in +% IEEE Transactions on Automatic Control, vol. 64, no. 7, pp. 2712-2726, July 2019, +% doi: 10.1109/TAC.2018.2867916. + import casadi.* lambda=SX.sym('lambda',nx,1); % the i th multiplier for equality constraints @@ -31,32 +40,52 @@ %% Generate some functions +% Computing the Jacobian of the implicit function [1, eqs.(1a)--(1e)]. A +% "lighther" description is also reported in [2, eqs.(1a)--(1e)]. For this +% part refers to the online documentation of CaSADI: https://web.casadi.org/docs/#document-function f_fun = Function('f_fun', {states,controls,params,alg}, {SX.zeros(nx,1)+x_dot},{'states','controls','params','alg'},{'xdot'}); +% Initializate the dimension of the symbolic vectors in CaSADI impl_jac_f_x = SX.zeros(nx,nx)+jacobian(impl_f,states); impl_jac_f_u = SX.zeros(nx,nu)+jacobian(impl_f,controls); impl_jac_f_xdot = SX.zeros(nx,nx)+jacobian(impl_f,xdot); +% Allocating the Functions in CaSADI. For this part refers to the online +% documentation: https://web.casadi.org/docs/#document-function impl_f_fun = Function('impl_f_fun',{states,controls,params,xdot,alg},{SX.zeros(nx,1) + impl_f}); impl_jac_f_x_fun = Function('impl_jac_f_x_fun',{states,controls,params,xdot,alg},{impl_jac_f_x}); impl_jac_f_u_fun = Function('impl_jac_f_u_fun',{states,controls,params,xdot,alg},{impl_jac_f_u}); impl_jac_f_xdot_fun = Function('impl_jac_f_xdot_fun',{states,controls,params,xdot,alg},{impl_jac_f_xdot}); +% In case the number of algebraic states is greater than zero if nz>0 + % F g_fun = Function('g_fun', {states,controls,params,xdot,alg}, {SX.zeros(nz,1)+z_fun},{'states','controls','params','xdot','alg'},{'zfun'}); + % Initializate the dimension of the symbolic vectors in CaSADI in the + % case the number of algebraic states is greater than zero impl_jac_f_z = SX.zeros(nx,nz)+jacobian(impl_f,alg); impl_jac_g_x = SX.zeros(nz,nx)+jacobian(z_fun,states); impl_jac_g_u = SX.zeros(nz,nu)+jacobian(z_fun,controls); impl_jac_g_z = SX.zeros(nz,nz)+jacobian(z_fun,alg); + % Allocating the Functions in CaSADI in the case the number of + % algebraic states is greater than zero. For this part refers to the + % online documentation of CaSADI: https://web.casadi.org/docs/#document-function impl_jac_f_z_fun = Function('impl_jac_f_z_fun',{states,controls,params,xdot,alg},{impl_jac_f_z}); impl_jac_g_x_fun = Function('impl_jac_g_x_fun',{states,controls,params,xdot,alg},{impl_jac_g_x}); impl_jac_g_u_fun = Function('impl_jac_g_u_fun',{states,controls,params,xdot,alg},{impl_jac_g_u}); impl_jac_g_z_fun = Function('impl_jac_g_z_fun',{states,controls,params,xdot,alg},{impl_jac_g_z}); else + % In the case the number of algrebraic states is equal to zero. In + % CaSADI the expression {0} refers to an empty function. For more + % information refers to the online documentation of CaSADI: + % https://web.casadi.org/docs/#document-function g_fun = Function('g_fun', {states,controls,params,xdot,alg}, {0},{'states','controls','params','xdot','alg'},{'zfun'}); + % Allocating the Functions in CaSADI in the case the number of + % algebraic states is equal to zero. For more info refers to the online + % documentation of CaSADI: https://web.casadi.org/docs/#document-function impl_jac_f_z_fun = Function('impl_jac_f_z_fun',{states,controls,params,xdot,alg},{0}); impl_jac_g_x_fun = Function('impl_jac_g_x_fun',{states,controls,params,xdot,alg},{0}); impl_jac_g_u_fun = Function('impl_jac_g_u_fun',{states,controls,params,xdot,alg},{0}); @@ -65,31 +94,65 @@ Sx = SX.sym('Sx',nx,nx); Su = SX.sym('Su',nx,nu); +% In the case the number of algebraic states is equal to zero if nz==0 + % This part implements the ADJ-RTI scheme (see [2, Sec.II]. For the + % "jacobian" function refers to + % https://web.casadi.org/docs/#calculus-algorithmic-differentiation. + % Instead, the "jtimes" function is described in the subsection 3.9.1. + % of the online documentation of CaSADI. + % + % jtimes = For calculating a Jacobian-times-vector product, the jtimes function + % – performing forward mode AD – is often more efficient than creating the full + % Jacobian and performing a matrix-vector multiplication: vdeX = SX.zeros(nx,nx); vdeX = vdeX + jtimes(x_dot,states,Sx); vdeU = SX.zeros(nx,nu) + jacobian(x_dot,controls); vdeU = vdeU + jtimes(x_dot,states,Su); vdeFun = Function('vdeFun',{states,controls,params,Sx,Su,alg},{vdeX,vdeU}); + % ERK stands for explicit Runge–Kutta integrator (see [1, Sec. V]) adjW = SX.zeros(nx+nu,1) + jtimes(x_dot, [states;controls], lambda, true); adj_ERK_fun = Function('adj_ERK_fun',{states,controls,params,lambda,alg},{adjW}); +% in the opposite case (i.e., the number of algrebraic states is greater +% than zero) else + % ERK stands for explicit Runge–Kutta integrator (see [1, Sec. V]) vdeFun = Function('vdeFun',{states,controls,params,Sx,Su,alg},{0,0}); adj_ERK_fun = Function('adj_ERK_fun',{states,controls,params,lambda,alg},{0}); end %% objective and constraints +% The inner function and the path constraint function are described in [2]. +% For more information regarding the meaning, please refer to the following +% papers: +% - [3] C.Kirches, L. Wirsching, H. G. Bock, J. P. Schlöder, "Efficient direct +% multiple shooting for nonlinear model predictive control on long +% horizons", Journal of Process Control, Volume 22, Issue 3, March 2012, +% Pages 540--550, doi: 10.1016/j.jprocont.2012.01.008 +% - [4] R. Quirynen, M. Vukov, M. Zanon, M. Diehl, "Autogenerating microsecond +% solvers for nonlinear MPC: A tutorial using ACADO integrators", Optimal +% Control, doi: https://doi.org/10.1002/oca.2152 +% - [5] E. Rossi, M. Bruschetta, R. Carli, Y. Chen and M. Farina, "Online Nonlinear +% Model Predictive Control for tethered UAVs to perform a safe and constrained maneuver," +% 2019 18th European Control Conference (ECC), 2019, pp. 3996-4001, doi: 10.23919/ECC.2019.8796032. +% - [6] T. Albin, D. Ritter, D. Abel, N. Liberda, R. Quirynen and M. Diehl, +% "Nonlinear MPC for a two-stage turbocharged gasoline engine airpath", 2015 +% 54th IEEE Conference on Decision and Control (CDC), Osaka, 2015, pp. 849-856. +% doi: 10.1109/CDC.2015.7402335 +% Some numerical examples are also presented in [1] h_fun=Function('h_fun', {states,controls,params}, {h},{'states','controls','params'},{'h'}); hN_fun=Function('hN_fun', {states,params}, {hN},{'states','params'},{'hN'}); path_con_fun=Function('path_con_fun', {states,controls,params}, {general_con},{'states','controls','params'},{'general_con'}); path_con_N_fun=Function('path_con_N_fun', {states,params}, {general_con_N},{'states','params'},{'general_con_N'}); +% This part refers to [2, eq.(4)]. gxi = jacobian(obji,states)' + SX.zeros(nx,1); gui = jacobian(obji,controls)' + SX.zeros(nu,1); gxN = jacobian(objN,states)' + SX.zeros(nx,1); +% Gauss-Newton Hessian approximation [2, eq.(5)]. Hz = hessian(obji_GGN, aux) + SX.zeros(ny, ny); HzN = hessian(objN_GGN, auxN) + SX.zeros(nyN, nyN); @@ -103,14 +166,17 @@ Cui = jacobian(general_con, controls) + SX.zeros(nc, nu); CxN = jacobian(general_con_N, states) + SX.zeros(ncN, nx); +% Objective functions at the terminal stage (objN_fun) and not) obji_fun = Function('obji_fun',{states,controls,params,refs,Q},{obji+SX.zeros(1,1)}); objN_fun = Function('objN_fun',{states,params,refN,QN},{objN+SX.zeros(1,1)}); +% Creating functions in CaSADI Hz_fun = Function('Hz_fun',{aux,params,refs,Q},{Hz}); HzN_fun = Function('HzN_fun',{auxN,params,refN,QN},{HzN}); Hw = Hz_fun(h,params,refs,Q) + SX.zeros(ny,ny); HwN = HzN_fun(hN,params,refN,QN)+ SX.zeros(nyN,nyN); +% Creating functions in CaSADI Hi_fun=Function('Hi_fun',{states,controls,params,refs,Q},{Hw}); HN_fun=Function('HN_fun',{states,params,refN,QN},{HwN}); diff --git a/README.md b/README.md index 761da9d..a4db9ca 100644 --- a/README.md +++ b/README.md @@ -44,3 +44,29 @@ Install Xcode from app store ### To use HPIPM as the QP solver please read /doc/HPIPM-tutorial for the detailed installation process + +### References + +``` +@ARTICLE{8451938, + author={Chen, Yutao and Bruschetta, Mattia and Cuccato, Davide and Beghi, Alessandro}, + journal={IEEE Transactions on Automatic Control}, + title={An Adaptive Partial Sensitivity Updating Scheme for Fast Nonlinear Model Predictive Control}, + year={2019}, + volume={64}, + number={7}, + pages={2712-2726}, + doi={10.1109/TAC.2018.2867916}} +``` + +``` +@INPROCEEDINGS{8795788, + author={Chen, Yutao and Bruschetta, Mattia and Picotti, Enrico and Beghi, Alessandro}, + booktitle={2019 18th European Control Conference (ECC)}, + title={MATMPC - A MATLAB Based Toolbox for Real-time Nonlinear Model Predictive Control}, + year={2019}, + volume={}, + number={}, + pages={3365-3370}, + doi={10.23919/ECC.2019.8795788}} +``` diff --git a/doc/.gitignore b/doc/.gitignore new file mode 100644 index 0000000..b04e459 --- /dev/null +++ b/doc/.gitignore @@ -0,0 +1,17 @@ +#### LaTeX #### + +*.pdf +!**/figures/*.pdf +*.aux +*.log +*.tex.bak +*.bak +*.gz +*.out +*.dvi +*.bst +*.blg +*.tdo +*.toc +documentation.pdf +documentation.bbl