Description

Exercises for Lecture Course on Numerical Optimal Control (NOC)
Exercise 6: Optimal Control Formulations
Prof. Dr. Moritz Diehl, Andrea Zanelli, Dimitris Kouzoupis, Florian Messerer

Figure 1: Sketch of the pendulum with its position defined by the angle θ.
In this exercise we will have a look at the special structures of the optimal control formulations discussed in the class. To this end, we will consider the control of a pendulum, as illustrated in Figure 1. The pendulum position is described by the angle θ with corresponding angle velocity w and we can control it by applying a torque τ. This system can be described by idealized dynamics
θ˙ = ω,
(1) ω˙ = sin(θ) + τ,
where for simplicity we are ignoring all units. Defining the state x = (θ,ω), with control u = τ, our control aim is to swing up the pendulum to the upright position and balancing it there, corresponding to a state of xref = (0,0). We start with the pendulum hanging down in a rest position, i.e., x¯ = (π,0). We can express this as the discrete time optimal control problem
(2)
s.t. x0 = ¯x0 xi+1 = F(xi,ui), i = 0,…,N − 1,
where x0,…,xN ∈R2, u0,…,uN−1 ∈R, and F(x,u) describes the discretized dynamics obtained by applying the explicit RK4 integration scheme to (1). The objective expresses our desire to bring the state x close to the origin, while using not too much control effort, as quantified by deviations of the control variable from 0. We will consider different reformulations of problem (2) and analyze their structure.
1. Sequential approach: The state variables in (2) can be eliminated by means of a forward simulation such that a more compact optimization problem is obtained,
min Φ(U), (3)
U
1
with U = (u0,…,uN−1). Using CasADi, formulate Φ(U) and set up problem (3), where F(xi,ui) is an RK4 integrator with 10 intermediate RK4 steps over a total integration time of h = 0.1. Use N = 50 for the discrete time horizon (such that the continuous time horizon is T = 5) and ¯x0 = (π, 0) as initial state. Solve the problem using IPOPT and plot the state and control trajectories in two separate plots. Use ui = 0.1, i = 0,…,N − 1, as initial guess. Before solving the problem, compute the Hessian of Φ(u) for some non-zero u and plot its sparsity pattern using the MATLAB command spy.
Hint: You can use rk4step.m provided with this exercise to do one step of the RK4 scheme.
2. Simultaneous approach: If both states and controls are kept as decision variables, a larger problem is obtained with a special sparsity structure that can be exploited by tailored algorithms.
(a) Using CasADi and the provided template simultaneous.m, formulate problem (2) without eliminating the state variables, solve it with IPOPT and plot the state and control trajectories. Make sure that you obtain the same result as with the sequential approach. Eliminate the initial state, by fixing it to x0 = ¯x0 and use ui = 0.1, i = 0,…,N − 1 and xi = (0.1,0.1), i = 1,…,N as initial guess. Before solving the problem compute the Hessian of the Lagrangian and plot its sparsity pattern using the MATLAB command spy. In order to highlight the block-banded structure of the Hessian, reorder the variables as [u0, λ1, x1, u1,λ2, …]. For this last part, you can set the equality multipliers to any arbitrary non-zero value since we are only interested in the sparsity pattern here.
(b) Implement now your own full-step Newton method and use it to solve problem (2). Initialize both primal and dual variables to 0.1. Terminate the algorithm when the 2norm of the right-hand-side is less than 10−6. Plot the resulting trajectories and make sure you obtain the same results as with the previous implementations.
Hint: you already have the Lagrangian of the problem from a) which can be used to compute the right-hand-side of the linear system associated with the Newton steps. Analogously, the matrix coefficient of the linear system can be obtained by computing the Hessian of the Lagrangian.
(c) Make a copy of your script from a) and modify the formulation in order to include box constraints on the inputs:
− 1 ≤ ui ≤ 1, ∀i = 0,…, N − 1. (4)
Solve the modified problem with IPOPT using the same initial guess as before and plot the obtained trajectories.
(d) [Bonus] Make a copy of your script from b) and modify the formulation by including logarithmic barrier terms in the cost that approximate the constraints introduced in c):
, (5)
where ρ = 0.01 is a fixed parameter. Solve the formulation with your implementation of the Newton method. In order to make sure that the arguments of the logarithms in the cost are always positive, implement a line search: starting from a full step, check whether any of the constraints are violated and half the step size until the arguments of the logarithms are positive. Initialize both primal and dual variables to 0.1.
2