Description

Exercise 3: Equality Constrained Optimization
Prof. Dr. Moritz Diehl, Andrea Zanelli, Dimitris Kouzoupis, Florian Messerer

In this sheet we will build on the previous exercise by implementing a Newton-type algorithm for equality constrained problems and looking into linear independence constraint qualification.
1. Newton method for equality constrained problems. Consider the following equality constrained optimization problem:
(1a) s.t. g(x,y) := x + (1 − y)2 = 0. (1b)
In this exercise we will implement a simple Newton-type algorithm that can be used to solve problem (1).
1 x = MX.sym(‘x’,2,1);
2 h = Function(‘h’, {x}, {x(1)ˆ2 + x(2)})
3 h = returntypes(‘full’, h);
(a) Compute on paper the gradients of f and g and their Hessian.

(b) Write on paper the Karush-Kuhn-Tucker (KKT) conditions for problem (1). Are these conditions necessary for optimality? Are they sufficient?
Lagrangian for this problem: L(x,y,λ) = f(x,y) + λ>g(x,y), λ ∈R KKT conditions (first order necessary conditions):
∇(x,y)L(x∗,y∗,λ∗) = ∇f(x∗,y∗) + ∇g(x∗,y∗)λ∗ = 0 g(x∗,y∗) = 0
If LICQ holds at (x∗,y∗), the KKT conditions are necessary conditions for optimality, but in general not sufficient. If the problem is additionally convex, they are also sufficient. The given specific problem is not convex, therefore they are not sufficient.
(c) In the provided template implement f and g and their Jacobians and Hessians as CasADi functions. With returntypes() you can set the output type to a normal MATLAB value. Like this you do not have to repeatedly use full(), which can be convenient sometimes.
(d) The KKT conditions derived in (b) can be written in compact form as
r(w) = 0, (2)
where w := [x,y,λ]T and λ is the Lagrange multiplier associated with the equality constraint g(x,y) = 0. Using the template provided, implement the following Newton-type method:
wk+1 = wk − M−1r(wk), (3)
where M ≈∇r(wk) is an approximation of the exact Jacobian of r. Test your implementation with two different Hessian approximations: i) B = ρI2 for different values of ρ and ii) B = ∇2f(xk,yk) + λ∇2g(xk,yk). Initialize the iterates at w0 = [1,−1,1]T and run the algorithm for N = 100 iterations. Plot the iterates in the x−y space. When using the fixed Hessian approximation, does the algorithm converge for ρ = 100? And for ρ = 600?
Divergence for ρ = 100, slow convergence for ρ = 600.
2. Linear independence constraint qualification. Consider the problem of finding the optimal way of throwing a ball such that progress in the horizontal coordinate after a fixed time T is maximized. The dynamics of the system can be modeled in two dimensions by the following differential equation: p˙y = vy,

where py and pz represent the y and z coordinate of the ball respectively and vy and vz the components of its velocity. The ball is subject to drag force with drag coefficient d, side wind w and gravitational acceleration g. In order to achieve the desired goal, we formulate the following optimization problem:
min
v − py(T;v) (4a)
s.t. −pz(T;v) ≤ 0, (4b)
−α(py(T;v) − 10) − pz(T;v) ≤ 0, (4c)
, (4d)
where v := [vy(0), vz(0)]T are the decision variables and p(T;v) is the output of an RK4 integrator that discretizes the dynamics of the system. Additional constraints have been added to the formulation that represent the requirement that the ball has to be above the ground at time T and that it has to be in the half-space defined by the linear constraint (4c), where α ∈ [−1, 1] is a fixed parameter.
(a) Implement the differential equation of the system in the provided ballistic dynamics.m.
(b) Using the provided template, implement and solve the optimization problem for different values of α ∈ [−1,1]. The template shows how the NLP constructed by CasADi can be parametrized by α. Like this the value of α can be conveniently changed without needing to construct a new NLP each time. Plot the normalized gradients of the constraints as three vectors. What happens when α = 0? For α ≥ α¯ for some ¯α ≥ 0 the problem becomes infeasible. What happens to the three vectors as α approaches ¯α? It is not required to compute the exact value of ¯α.
at α = 0: constraints (4b) and (4c) are identical. Therefore their gradients are parallel everywhere, which leads to LICQ violation.
at α = ¯α: Exactly one feasible point is left, which is therefore also the optimal value. At this point constraint (4c) and (4d) are tangential to each other, i.e., their gradients are parallel. Therefore LICQ is violated.