A simple example of sensitivity analysis

In this example we explore the parametric sensitivity of an NLP problem. We will use the csnlp.wrappers.NlpSensitivity wrapper to achieve so, but we’ll skip over some of the formalities of the library, as we assume you are already a bit familiar with it. If this is not the case, for more introductory examples on csnlp, see the other Introductory examples.

The problem is inspired the example from this CasADi blog post. The following problem is considered

\[\begin{split}\begin{aligned} \min_{x} \quad & (1 - x_1)^2 + p_1 (x_2 - x_1^2)^2 \\ \text{s.t.} \quad & x_1 \ge 0 \\ & \frac{p_2^2}{4} \le (x_1 + 0.5)^2 + x_2^2 \le p_2^2. \end{aligned}\end{split}\]

where \(p = [p_1, p_2]^\top\) are the parameters of the problem.

We start with the usual imports.

import casadi as cs
import matplotlib.pyplot as plt
import numpy as np

from csnlp import Nlp, wrappers

Building the NLP

We then build the NLP in the usual way. We create the variables and the parameters of the problem, and then we add the objective function and the constraints. Lastly, we initialize the solver.

nlp = Nlp[cs.MX](sym_type="MX")
x = nlp.variable("x", (2, 1), lb=[[0], [-np.inf]])[0]
p = nlp.parameter("p", (2, 1))

nlp.minimize((1 - x[0]) ** 2 + p[0] * (x[1] - x[0] ** 2) ** 2)
g = (x[0] + 0.5) ** 2 + x[1] ** 2
nlp.constraint("c1", (p[1] / 2) ** 2, "<=", g)
_, lam = nlp.constraint("c2", g, "<=", p[1] ** 2)
nlp.init_solver({"print_time": False, "ipopt": {"sb": "yes", "print_level": 0}})

To spice things up, we’ll convert the NLP to a casadi.Function via the method csnlp.Nlp.to_function. This streamlines the process of evaluating the NLP multiple times. The function takes as input the parametrization and the initial guess for the primal variables, and returns the optimal primal solution.

M = nlp.to_function("M", [p, x], [x], ["p", "x0"], ["x"])

We can call M to solve the optimization for different values of \(p\).

p_values = [[0.2, 1.2], [0.2, 1.45], [0.2, 1.9]]
X_opt = [M(p, 0.0).full() for p in p_values]

A nice visualization follows. In red we have the infeasible region, while the star marks the optimal solution. The contours represent the objective function’s levels.

_, axs = plt.subplots(1, 3, sharey=True, constrained_layout=True)

ts = np.linspace(0, 2 * np.pi, 100)
cos_ts, sin_ts = np.cos(ts), np.sin(ts)
X, Y = np.meshgrid(*[np.linspace(-1.5, 2, 100)] * 2)

for ax, (p1, p2), (x1_opt, x2_opt) in zip(axs, p_values, X_opt):
    ax.contour(X, Y, (1 - X) ** 2 + p1 * (Y - X**2) ** 2, 100, alpha=0.5)
    ax.fill(
        p2 / 2 * cos_ts - 0.5, p2 / 2 * sin_ts, facecolor="r", alpha=0.7, edgecolor="k"
    )
    ax.fill(
        np.concatenate([p2 * np.cos(ts - np.pi) - 0.5, 10 * np.cos(np.pi - ts) - 0.5]),
        np.concatenate([p2 * np.sin(ts - np.pi), 10 * sin_ts]),
        facecolor="r",
        alpha=0.7,
        edgecolor="k",
    )
    ax.fill(
        [-1.5, 0, 0, -1.5], [-1.5, -1.5, 2, 2], facecolor="r", alpha=0.7, edgecolor="k"
    )
    ax.plot(x1_opt, x2_opt, "*", markersize=12, linewidth=2)
    ax.set_aspect("equal", "box")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_xlim(-1.5, 2)
    ax.set_ylim(-1.5, 2)
    ax.set_title(rf"$p_1={p1:.2f}, \ p_2={p2:.2f}$")

plt.show()

From this plot, we can appreciate how the optimal solution changes as we vary the parametrization of the problem, especially when the constraint becomes inactive. This begs the question: how does the optimal solution vary along \(p_2\)? We can answer this question via sensitivity analysis.

Sensitivity analysis

To be as general as possible, consider a generic function that maps the optimal primal and dual variables and parameters to a vector of expressions. We call this mapping \(Z(x(p), \lambda(p), p)\). Note that \(x(p)\) and \(\lambda(p)\) are the functions of the parameters \(p\), as we just saw, but this relationship is implicit and most of the time difficult to express analytically. Our goal is to numerically compute its sensitivity w.r.t. \(p_2\), i.e.,

\[\frac{\partial Z}{\partial p_2} \quad \text{and} \quad \frac{\partial^2 Z}{\partial p_2^2}.\]

In order to to do, let us build a very generic function Zfcn (note that its entries have no clear meaning, are just arbritary continuous functions). Again, we can use the csnlp.Nlp.to_function method to convert it to a CasADi function for ease of evaluations.

z1 = lambda x: x[1, :] - x[0, :]
z2 = lambda x, lam, p: (x[1, :] ** p[1] - x[0, :]) * cs.exp(-10 * lam + p[0]) / p[1]
z3 = lambda x, p: x[1, :] ** (1 / p[1]) - x[0, :]
z4 = lambda lam, p: cs.exp(-10 * lam + p[0]) / p[1]
Z = cs.vertcat(
    z1(x),
    z2(x, lam, p) ** 2,
    1 / (1 + z3(x, p)),
    z4(lam, p),
    z3(x, p) * (-1 - 10 * z2(x, lam, p)),
    z4(lam, p) / (1 + z1(x)),
)
Zfcn = nlp.to_function("Z", [p, x], [Z], ["p", "x0"], ["z"])

To enable the computation of the parametric sensitivities, we need to wrap the NLP in a csnlp.wrappers.NlpSensitivity object. This takes in the target parameter \(p_2\), because that’s the parameter we are interested in.

nlp_ = wrappers.NlpSensitivity[cs.MX](nlp, target_parameters=p[1])

The symbolic representation of the parametric sensitivities (jacobian and hessian) is obtained via the csnlp.wrappers.NlpSensitivity.parametric_sensitivity method as follows. By default, the method does not compute second order information, because it can be very costly, and, if no expression is passed, will compute the sensitivity of the primal variables w.r.t. the target_parameters.

J, H = nlp_.parametric_sensitivity(expr=Z, second_order=True)

Once the symbols are available, we can evaluate them numerically after a solution is computed. We do so for different values of \(p_2\), and store the results in some lists for later plotting.

Z_values, J_values, H_values = [], [], []
for p_value in p_values:
    sol = nlp.solve(pars={"p": p_value})
    Z_values.append(sol.value(Z).full().flatten())  # also evaluate the actual function
    J_values.append(sol.value(J).full().flatten())
    H_values.append(sol.value(H).full().flatten())

Visualization of the sensitivities

Before jumping into the visualization of the sensitivities themselves, let’s compute the function \(Z\) for a range of realizations of \(p_2\) (with \(p_1=0.2\)). This will also serve as a baseline to understand whether the sensitivities follow the actual plot.

N = 300
p_values_all = np.vstack((np.full(N, 0.2), np.linspace(1, 2, N)))
z_values_all = np.concatenate([Zfcn(p, 0.0) for p in p_values_all.T], axis=1)

We can now plot the function \(Z\) for the different values of \(p_2\). After that, we can plot the parametric sensitivities. To do so, we plot the second order Taylor expansion that approximates the function \(Z\) around the optimal solution for different values of \(p_2\). This approximation obviously makes use of the gradient and hessian information of the function at that point.

fig, axs = plt.subplots(3, 2, sharex=True, constrained_layout=True)
for i, ax in enumerate(axs.flat):
    ax.plot(p_values_all[1], z_values_all[i], "k-", lw=4)
    ax.set_ylabel(f"$Z_{i}$")
    if i in (4, 5):
        ax.set_xlabel("$p_2$")
    ax.set_ylim(ax.get_ylim() + np.asarray([-0.1, 0.1]))

t = np.linspace(1, 2, 100)
for i in range(len(p_values)):
    p = p_values[i][1]
    Z_value = Z_values[i]
    J_value = J_values[i]
    H_value = H_values[i]
    for ax, z, j, h in zip(axs.flat, Z_value, J_value, H_value):
        ax.plot(p, z, "o", color=f"C{i}", markersize=6)
        ax.plot(t, z + j * (t - p) + 0.5 * h * (t - p) ** 2, color=f"C{i}", ls="--")

plt.show()