Skip to content

Commit

Permalink
Added tutorial for the non-sparse solvers
Browse files Browse the repository at this point in the history
  • Loading branch information
@pablokrupa
pablokrupa committed Sep 3, 2021
1 parent be0ea82 commit 5f15f5f
Show file tree
Hide file tree
Showing 2 changed files with 183 additions and 11 deletions.
23 changes: 12 additions & 11 deletions examples/basic_tutorial.m
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,32 @@
%% Welcome to SPCIES: Suite of Predictive Controllers for Industrial Embedded systems.
%
%
% This is a basic introductory tutorial to the main use of SPCIES: the generation
% of tailored solvers for different MPC formulations.
%
%
% The toolbox can create solvers in various programming languages.
% This tutorial will guide you through the creation of a MEX function for Matlab
% for a standard MPC formulation that we label 'laxMPC', which is given by:
%
%
% min_{x, u} sum_{i = 0}^{N-1} \| x_i - xr \|^2_Q + \| u_i - ur \|^2_R + \| x_N - xr \|^2_T
%
% s.t. x_0 = x(k)
% x_{i+1} = A x_i + B u_i, i = 0...N-1
% LBx <= x_i <= UBx, i = 1...N-1
% LBu <= u_i <= UBx, i = 0...N-1
%
% LBu <= u_i <= UBu, i = 0...N-1
%
% where x_i and u_i are the predicted states and control actions, respectively;
% x_r and u_r are the state and input reference, respectively; Q, R and T are
% the positive definite cost function matrices; LBx, LBu are the lower bounds
% for the state and input, respectively; UBx, UBu the upper bounds; and
% N is the prediction horizon.
%
% For additional information about the formulation, we refer the user to:
%
%
% P. Krupa, D. Limon, T. Alamo, "Implementation of model predictive control in
% programmable logic controllers", Transactions on Control Systems Technology, 2020.
%
% Specifically, this formulation is given in equation (9) of the above reference.
%
%
clear; clc;

%% STEP 1: Problem definition.
Expand Down Expand Up @@ -73,14 +73,14 @@
% Once again, the name of the fields must be as shown.
param = struct('Q', Q, 'R', R, 'T', T, 'N', N);

%% STEP 3: Generate the solver
%% STEP 3: Generate the solver.
% In this step, we will call the main function of SPCIES to generate a
% MEX file containing the sparse solver.
% In this tutorial we will create an ADMM-based solver for the laxMPC formulation.

% Before we generate the solver, we must sets its options. This is not strictly
% necessary, since default values are provided, but it is recommended. Especially
% for the penalty parameter 'rho' of ADMM, since performance is highly dependant
% for the penalty parameter 'rho' of ADMM, since performance is highly dependent
% on its value.
% To set the options we must create a structure with the appropriate fields.
% We will set some options. For a full list of the available options, please
Expand All @@ -103,7 +103,7 @@
% To save the files into the current working directory use options.directory = './';

% Since we are saving the solver into the default directory $SPCIES$/generated_solvers,
% we recomment that you clear the directory of all previous controllers.
% we recommend that you clear the directory of all previous controllers.
% Controllers are overwritten by default, but it is best to be safe that sorry.
% You can clear $SPCIES$/generated_solvers by running:
spcies_clear;
Expand All @@ -130,7 +130,7 @@
% If using the Linux version of Matlab, SPCIES currently assumes that the
% gcc compiler is being used.

%% STEP 4: Use the MEX file
%% STEP 4: Use the MEX file.
% To show how to use the generated MEX file, we will perform a closed-loop test.

% First, we set the conditions of the test
Expand Down Expand Up @@ -207,3 +207,4 @@
xlabel('Sample time');
ylabel('Number of iterations');
grid on;

171 changes: 171 additions & 0 deletions examples/tutorial_non_sparse.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
%% Tutorial: The non-sparse Matlab solvers
%
% This tutorial covers the basics of the non-sparse version of the solvers
% available in the SPCIES toolbox.
%
% Before reading through this tutorial we recommend you first take a look
% at the basic tutorial (basic_tutorial.m).
%
% The main objective of the SPCIES toolbox y to generate sparse tailored
% solvers for different MPC formulations. These solvers are based on
% optimized code that is difficult to understand, to change or to extract
% additional information from (beyond the outputs of the function, that is).
%
% To alleviate this issue, the SPCIES toolbox includes non-sparse versions of
% all its solvers written in plan Matlab code. These solvers, which are
% stored in the folder $SPCIES$/platforms/Matlab, where $SPCIES$ is the root
% directory of the toolbox (call spcies_get_root_directory), provide several
% advantages and disadvantages with respect to the sparse solvers generated
% using the spcies_gen_controller function:
%
% ADVANTAGES:
% - The non-sparse solvers return more information than the sparse versions.
% For instance, they can be asked to return the historic of the values
% of the decision variables during its iterations (see below for how
% to do this).
% - The Matlab code is easier to follow and understand than the sparse code.
% This allows the user to follow along with what the algorithm is doing,
% especially if it is compared to the articles and other documentation
% explaining the solvers. It also helps the user to understand what the
% sparse code is doing.
% This is only useful if your are interested in understanding the solver.
% - It facilitates making changes to the algorithm to test new ideas or suit
% your particular needs. This is interesting if you are doing research or
% need access to data that is hard to extract from the sparse solver.
% - It allows for debugging using Matlab's debugging tools.
% - It does not require a compiler to generate the solver. It directly runs
% in Matlab.
%
% DISADVANTAGES:
% - The non-sparse solver is significantly slower that the sparse version.
% Especially so if repeatedly called, since it performs all the offline
% computations each times it is called, whereas the sparse version only
% performs them when the spcies_gen_controller function is called.
% - It can only be executed in Matlab, whereas the sparse solvers can be
% generated for different programming languages and target platforms.
%
% In short, these functions are useful if you want to be able to debug, make
% changes to the algorithm or need data not returned by the sparse version.
%
% This tutorial will show an example of how to use one of these functions.
% The usage of the others is very similar.
% The reader will not that the procedure is very similar to the one used for
% the sparse versions (see basic_tutorial.m).
%
% In this tutorial we will use the solver spcies_MPCT_EADMM_solver, which
% is a solver, based on the EADMM algorithm, for the MPC formulation:
%
% min_{x, u, xs, us} sum_{i=0}^{N-1} \| x_i - xs \|^2_Q + \| u_i - us \|^2_R +
% + \| xs - xr \|^2_T + \| us - ur \\^2_S
%
% s.t. x_0 = x(k)
% x_{i+1} = A x_i + B u_i, i = 0...N-1
% LBx <= x_i <= UBx, i = 1...N-1
% LBu <= u_i <= UBu, i = 0...N-1
% xs = A xs + B us
% LBx + epsilon_x <= xs <= UBx - epsilon_x
% LBu + epsilon_u <= us <= UBu - epsilon_u
% x_N = xs
%
% where x_i and u_i are the predicted states and control actions, respectively;
% x_r and u_r are the state and input reference, respectively; xs and us are
% artificial reference, Q, R, T and S, the positive definite cost function matrices;
% LBx, LBu are the lower bounds for the state and input, respectively; UBx, UBu the
% upper bounds; epsilon_x and epsilon_u vectors with arbitrarily small components
% and N is the prediction horizon.
%
% For additional information about the formulation, which we label MPCT, we refer
% the reader to (and the references therein):
%
% "Implementation of model predictive control for tracking in embedded systems
% using a sparse extended ADMM algorithm", by P. Krupa, I. Alvarado, D. Limon
% and T. Alamo, arXiv preprint: 2008:09071v2, 2020.
%
clear; clc;

%% STEP 1: Problem definition.
% We consider the same system of connected masses that we used in the basic_tutorial.

p = 3; % Number of objects
M = [1; 0.5; 1]; % Mass of each object
K = 2*ones(p+1, 1); % Spring constant of each spring
F = [1; zeros(p-2, 1); 1]; % Objects in which an external force is applied

sysC = gen_oscillating_masses(M, K, F); % Generate the continuous-time model

n = size(sysC.A, 1); % For convenience, we save the size of the state dimension
m = size(sysC.B, 2); % and of the input dimension

% We obtain a discrete-time state space model
Ts = 0.2; % Select the sample time
sysD = c2d(sysC, Ts); % Discrete-time state space model

% We set some bounds for the state and input
LBx = -[ones(p, 1); 1000*ones(p, 1)]; % Lower bound for the state
UBx = [0.3; 0.3; 0.3; 1000*ones(p, 1)]; % Upper bound for the state
LBu = -[0.8; 0.8]; % Lower bound for the input
UBu = [0.8; 0.8]; % Upper bound for the input

% The ingredients that define the state space model have to be saved into a structure.
% The name of the fields of the structure have to be the ones shown here.

sys = struct('A', sysD.A, 'B', sysD.B, 'LBx', LBx, 'UBx', UBx, 'LBu', LBu, 'UBu', UBu);

%% STEP 2: Design the MPC controller.
% We need to determine the values of the ingredients of the MPC controller.

Q = blkdiag(15*eye(p), 1*eye(p)); % Remember that Q and R must be positive definite
R = 0.1*eye(m);
T = 10*Q; % In this formulation is it typically beneficial to make T larger than Q
S = R;
N = 10;

% The ingredients of the controller have to be saved into a structure.
% Once again, the name of the fields must be as shown.
param = struct('Q', Q, 'R', R, 'T', T, 'S', S, 'N', N);


%% STEP 3: Select options for the solver (optional).

% Set the options for the solver
options.rho_base = 2; % Base penalty parameter of the EADMM algorithm
options.rho_mult = 20; % Multiplication factor for the base penalty parameter
options.k_max = 5000; % Maximum number of iterations of the solver
options.tol = 1e-3; % Exit tolerance of the solver

% Notice that in this case there are no options for the toolbox. Only for the solver.

%% NOTICE: Steps 1 and 2 are exactly the same as for the sparse solver.
% The selection of the options is also the same.
% The only difference is that you do not have to generate the
% solver, you only need to call it, as you will now see.

%% STEP 4: Use the solver.
% We now show how to use the solver. We will only show a single call, but
% its use for a closed-loop simulation follows similarly.

% We need to determine the current state of the system and the reference.
x = zeros(n, 1); % Current state of the system
ur = 0.5*ones(m, 1); % We select a reference ur
xr = (sys.A - eye(n))\(-sys.B*ur); % We compute a state reference xr that is a steady state

% We call the solver
[u, k, e_flag, Hist] = spcies_MPCT_EADMM_solver(x, xr, ur, 'sys', sys, 'param', param,...
'options', options, 'genHist', 2);

% All the non-sparse solvers require the three inputs (x, xr, ur), the 'sys' structure
% and the 'param' structure. The 'options' structure is optional.
% They also accept other additional arguments. In particular, an important one to know
% about is 'genHist'. This argument can take the values 0, 1 or 2. Its value determines
% the amount of information that is saved and returned inside the structure 'Hist'.
% Lower values of 'genHist' return less information, but may have shorter execution times.

% The outputs of the non-sparse solvers are always the same:
% - u: Control input to apply to the system.
% - k: Number of iterations of the algorithm.
% - e_flag: Exit flag of the algorithm. Positive integers indicate success and negative
% numbers indicate that some problem has occurred. For the exact flags and their
% meaning we refer the reader to the help of each solver.
% - Hist: structure containing information about the solution and iterates of the solver.
% The extent of the information gathered is determined by the input 'genHist'.
% For an exact list of its fields we refer the reader to the help of each solver.

0 comments on commit 5f15f5f

Please sign in to comment.