Description

Exercise 5: Algorithmic Differentiation
Prof. Dr. Moritz Diehl, Dimitris Kouzoupis, Andrea Zanelli, Florian Messerer

The aim of this exercise is to gain experience with the two modes of algorithmic differentiation (AD) discussed in the class.
1. Forward and backward algorithmic differentiation: Consider the following discrete-time dynamical system:
xk+1 = xk + h((1 − xk)xk + uk), (1)
with state xk ∈R, input uk ∈R and h ∈R+ a constant parameter (you can think of it as the time step of an explicit Euler integrator). We are interested in simulating the dynamics forward for N steps, starting from the initial value x0 = ¯x0 and computing the derivatives of the obtained states with respect to the controls:
(2)
(a) Fix ¯x0 = 0.5, N = 50, h = 0.1. Make sure to define them once only, so you can easily adapt their values later. Using CasADi, implement the function Φ : RN →RN that maps controls to the obtained state trajectory
x = Φ(u), (3)
where x = (x1,…,xN) and u = (u0,…,uN−1) denote the vector of stacked states and controls respectively. Define a CasADi function that outputs the Jacobian of x with respect to u,
. (4)
You will use the output of this function as a reference for your implementations in the rest of the exercise.
(b) Implement a MATLAB function forw AD(u, m, x0, h) that takes as input a vector containing the values for u and a scalar m and returns the derivative , i.e., the m-th column of the Jacobian, evaluated at input u, using forward AD. Check that the result provided by your implementation is equal to the corresponding entries in the output obtained with CasADi, e.g., by evaluating both at randomly generated values of u, u tst = rand(N,1), and comparing the maximum of their elementwise absolute difference. You should be able to reach machine precision, i.e., order of magnitude around 10−16.
(c) Analogously, implement a MATLAB function back AD that takes as input u, a scalar m as well as the parameters. It returns the derivative , i.e., the m-th row of the Jacobian, using backward AD. Check that the result provided by your implementation is equal to the corresponding entry in the output obtained with CasADi.
(d) Implement now a function J FAD that takes as inputs u and a scalar m and, using forward
AD, computes the last m rows of the Jacobian containing the derivatives of the last m states in the simulation with respect to the all the controls. Again, validate your results against the reference output.
Note: You can just call forw AD several times, but then you unnecessarily repeat some computations. Can you do better?
1
(e) Analogously, implement a function J BAD that takes as inputs u and a scalar m and computes the last m rows of the Jacobian using your implementation of backward AD and validate it against the reference.
Again: You can just call back AD several times, but then you would unnecessarily repeat some computations. Can you do better?
(f) Which of the two implementation do you expect to be more performant for small values of m? Which one for high values of m? Why?
(g) Run your implementations for m ranging from 1 to N and measure the execution time using the MATLAB functions tic and toc. For this simulation choose h = 0.01 and N = 500. Plot the obtained execution times as a function of m. Do the results validate your considerations from the previous question?
2