%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (CHEM-E7225) ADVANCED PROCESS CONTROL D % #00 - Introduction to CasADi (Part 1) % % CasADi is a tool to quickly implement and efficiently solve dynamic % optimisation (optimal control) problems. You can also think of it as % a general toolkit for gradient-based optimisation % % Casadi is available in C++, Python, and Matlab/Octave % (There used to be an interface to Haskell, too) % % The numerical back-end for optimisation is IPOPT, plus others % The numerical back-end for simulation is SUNDIALS % % CasADI is open-source, LGPL-licensed % Download from https://casadi.org (current version: 3.5.5) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all; clear; clc addpath('~/Documents/MATLAB/casadi') % --> addpath('~/change_me_with_path_to/casadi-v3.5.5') import casadi.* % CasADi exists because derivatives are central in gradient-based optimisation, % yet they are extremely tedious and time-consuming to compute % % 1. By hand % 2. Symbolic diffentiation % 3. Finite, numerical, differences % 4. >> Algorithmic differentiation, AD << % CasADi has a symbolic framework for defining symbolic expressions using the everything-is-matrix % Matlab-alike style. All matrixes are defined to be sparse for storage efficiency %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SCALAR EXPRESSIONS (SX) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create two scalar symbolic primitives called "x" and "y" % % > Different variables can have the same name % > The identifier is the return x = SX.sym('x'); % A 1x1 matrix containing the symbol called "x" y = SX.sym('y'); % A 1x1 matrix containing the symbol called "y" whos x y %% Combine the symbols into an simple expression using common math operators % Visualise the expression using 'disp' or 'display' z = x * sin(x + y); % An expression z: R^2 -> R display(z) whos z %% Expressions can also be formed by composition of several (sub)expressions % For instance, consider the expression z = 10*exp((x + 10y^2) * sqrt(y)) z1 = x + 10*y^2; % An expression z_1: R^2 -> R, z_1(x,y) = x + 10y^2 z2 = z1 * sqrt(y); % An expression z_2: R^2 -> R, z_2(x,y) = z_1(x,y) * sqrt(y) = (x + 10y^2) * sqrt(y) z = 10 + exp(z2); % An expression z: R^2 -> R, z(x,y) = 10*exp(z_2(x,y)) = 10*exp((x + 10y^2) * sqrt(y)) display(z) whos z %% Expressions can also be formed by composition of several (sub)expressions >>iteratively<< % For instance, consider the expression z = 1 - x^2/2! + x^4/4! - x^6/6! + ... z = 0; for i = 0:2:8 z = z + (-1)^i * x^i / factorial(i); end display(z) whos z %% Compute the Jacobian of function z = x * sin(x + y) with respect to x % and with respect to y, then name them as Jx and Jy respectively z = x * sin(x + y); % An expression z: R^2 -> R Jz_x = jacobian(z, x); % dz(x,y)/dx, size (1 x 1) Jz_y = jacobian(z, y); % dz(x,y)/dy, size (1 x 1) display(Jz_x); display(Jz_y); whos Jz_x Jz_y % We can combine them to compute the Jacobian matrix, or use CasADi to do it tmpJz = [Jz_x Jz_y]; % [dz/dx dz/dy], size (1 x 2) Jz = jacobian(z, [x;y]); % [dz/dx dz/dy], size (1 x 2) display(tmpJz) display(Jz) whos tmpJz Jz % ... and we can also ask CasADi to compute the Gradient of z (the Jacobian transpose) Gz = gradient(z, [x;y]); % [dz/dx dz/dy]', size (2 x 1) display(Gz) whos Jz Gz %% Compute the Hessian of function z = x * sin(x + y) z = x * sin(x + y); % An expression z: R^2 -> R Hz = hessian(z, [x;y]); % |d^2z/dx^2 d^2z/dxdy| display(Hz) whos Hz %% Expressions can be combined into vector-valued expressions Z = [z, z^2, (z^2)^(1/2)]'; % A vector-valued expression Z: R^2 -> R^3 display(Z) whos Z % We can then use CasADi to compute Jacobians and Hessians JZ = jacobian(Z, [x;y]); % | dz_1/dx dz_1/dy| % | dz_2/dx dz_2/dy| % | dz_3/dx dz_3/dy| display(JZ) whos JZ HZ = hessian(Z, [x;y]); % This doesn't work! Why not? for i = 1:3 % In case needed, this would be HZ{i} = hessian(Z(i), [x;y]); % one possible way to compute Hessians end % for vector-valued expressions display(HZ) display(HZ{1}) whos HZ %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MATRIX EXPRESSIONS (SX and MX) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Create two `matrices' of scalar primitives called "A" and "b" and query % some of their properties using an appropriate member function and do % some other simple linear algebra stuff with them A = SX.sym('A',3,3); b = SX.sym('b',3); display(A) display(b) %% Matrix operations / Linear algebra normA = norm(A, 'fro'); % Frobenius-norm of matrix A, norm(A) = sum_{i=1}^N sum_{j=1}^N (A_ij)^2 traceA = trace(A); % Trace of matrix A, Tr(A) = sum_{i=1}^N A_ii A_x_b = A*b; % Matrix-vector multiplication hatx = solve(A,b); % Solves the linear system Ax = b (equivalent to A\b) %% Get and Set / Slicing A(1,1) = 0; % Set element A(1,1) A(2,2) = z; % Set element A(2,2) A(3,3) = SX.zeros(1,1); % Set element A(3,3) display(A) A_col = A(:,1); % Get the first column of matrix A A_row = A(1,:); % Get first row of matrix A display(A_col) display(A_row) %% CasADi provide a more efficient object for dealing with expressions with % matrix-valued operations, the Matrix eXpression (MX) type. % Consider the following operation using a matrix of SX objects: n = 3; A = SX.sym('A',n,n); B = SX.sym('B',n,n); C = A*B; % The operation is expanded into scalar operations % new expressions are created for each entry display(C) whos C % The same operation is represented in much more compact form using MX: A = MX.sym('A',n,n); B = MX.sym('B',n,n); C = A*B; % The operation is now represented by the % matrix multiplication explicitly display(C) whos C %% Matrix operations / Linear algebra A = MX.sym('A',n,n); b = MX.sym('b',n); normA = norm(A, 'fro'); % Frobenius-norm of matrix A, norm(A) = sum_{i=1}^N sum_{j=1}^N (A_ij)^2 traceA = trace(A); % Trace of matrix A, Tr(A) = sum_{i=1}^N A_ii A_x_b = A*b; % Matrix-vector multiplication hatx = solve(A,b); % Solves the linear system Ax = b (equivalent to A\b) %% Matrix-valued expressions are formed as in the scalar case (but with matrix operations!!) A = MX.sym('A',3,3); b = MX.sym('b',3); C = A * sin(b); % A (matrix) expression C: R^{3x3} x R^3 -> R^3 display(C) whos C %% Matrix-valued expressions can also be formed by composition A2 = A^2; % A (matrix) expression A_2: R^{3x3} x R^3 -> R^3, A_2(A,b) = A^2 b2 = sqrt(A2*b); % A (matrix) expression b_2: R^{3x3} x R^3 -> R^3, b_2(A,b) = sqrt(A2(A,b)*b) = sqrt(A^2 * b) C = A2\b2; % A (matrix) expression C: R^{3x3} x R^3 -> R^3, Solution of A2*x = b2 = (A^2)*x = sqrt(A^2 * b) display(C) whos C %% Expressions can also be formed by composition of several (sub)expressions >>iteratively<< % In this example, function C: R^{3x3} x R^{3x10} -> R given by C = \sum_{i=1}^10 norm(A * x_i)^2 A = MX.sym('A',3,3); x = MX.sym('x',3,10); C = 0; for i = 1:10 C = C + norm(A * x(:,i))^2; end display(C) whos C %% Compute the Jacobian/Gradient of a function z = x'*A*u with respect % to vectors x (10x1) and y (10x1), given matrix A (10x10). A = MX.sym('A',10,10); x = MX.sym('x',10); y = MX.sym('y',10); z = x'*A*y; % An expression z: R^(10x10) x R^10 x R^10 -> R Jz_xy = gradient(z, [x;y]); % Dz(A,x,y), size (? x ?) Jz_A = gradient(z, A); % dz(A,x,y)/dA, size (? x ?) display(Jz_xy) display(Jz_A) whos Jz_xy Jz_A % PS.: It is usally the case that compositions of matrix expressions can be % written explictly as one single expression in MATLAB, which can lead to % some improvemet in computational performance. For instance, which of the % above expressions can you write using a single line? %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % HANDS-ON EXERCISE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Write the following expresions using the MX type, % % x(1) = x0 % x(2) = A*x(1) + B*u(1) % x(3) = A*x(2) + B*u(2) % ... % x(N) = A*x(N-1) + B*u(N-1) % % given a 'horizon' N, matrices (A,B), and a symbol x0. (Play with the sizes!) % Then compute and inspect the Jacobian Jx = (dx / dx0) N = 10; nx = 2; nu = 1; A = MX.sym('A',nx,nx); B = MX.sym('B',nx,nu); u = MX.sym('u',nu,N-1); x0 = MX.sym('x0',nx); % [...] Complete the rest