.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/sensitivity/simple_sensitivity.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_sensitivity_simple_sensitivity.py: .. _simple_sensitivity_example: A simple example of sensitivity analysis ======================================== In this example we explore the parametric sensitivity of an NLP problem. We will use the :mod:`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 :mod:`csnlp`, see the other :ref:`introductory_examples`. The problem is inspired the example from `this CasADi blog post `_. The following problem is considered .. math:: \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} where :math:`p = [p_1, p_2]^\top` are the parameters of the problem. .. GENERATED FROM PYTHON SOURCE LINES 28-29 We start with the usual imports. .. GENERATED FROM PYTHON SOURCE LINES 29-36 .. code-block:: Python import casadi as cs import matplotlib.pyplot as plt import numpy as np from csnlp import Nlp, wrappers .. GENERATED FROM PYTHON SOURCE LINES 37-42 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. .. GENERATED FROM PYTHON SOURCE LINES 42-53 .. code-block:: Python 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}}) .. GENERATED FROM PYTHON SOURCE LINES 54-58 To spice things up, we'll convert the NLP to a :class:`casadi.Function` via the method :meth:`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. .. GENERATED FROM PYTHON SOURCE LINES 58-61 .. code-block:: Python M = nlp.to_function("M", [p, x], [x], ["p", "x0"], ["x"]) .. GENERATED FROM PYTHON SOURCE LINES 62-63 We can call ``M`` to solve the optimization for different values of :math:`p`. .. GENERATED FROM PYTHON SOURCE LINES 63-67 .. code-block:: Python 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] .. GENERATED FROM PYTHON SOURCE LINES 68-70 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. .. GENERATED FROM PYTHON SOURCE LINES 70-102 .. code-block:: Python _, 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() .. image-sg:: /auto_examples/sensitivity/images/sphx_glr_simple_sensitivity_001.png :alt: $p_1=0.20, \ p_2=1.20$, $p_1=0.20, \ p_2=1.45$, $p_1=0.20, \ p_2=1.90$ :srcset: /auto_examples/sensitivity/images/sphx_glr_simple_sensitivity_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 103-107 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 :math:`p_2`? We can answer this question via sensitivity analysis. .. GENERATED FROM PYTHON SOURCE LINES 109-126 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 :math:`Z(x(p), \lambda(p), p)`. Note that :math:`x(p)` and :math:`\lambda(p)` are the functions of the parameters :math:`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. :math:`p_2`, i.e., .. math:: \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 :meth:`csnlp.Nlp.to_function` method to convert it to a CasADi function for ease of evaluations. .. GENERATED FROM PYTHON SOURCE LINES 126-141 .. code-block:: Python 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"]) .. GENERATED FROM PYTHON SOURCE LINES 142-145 To enable the computation of the parametric sensitivities, we need to wrap the NLP in a :class:`csnlp.wrappers.NlpSensitivity` object. This takes in the target parameter :math:`p_2`, because that's the parameter we are interested in. .. GENERATED FROM PYTHON SOURCE LINES 145-148 .. code-block:: Python nlp_ = wrappers.NlpSensitivity[cs.MX](nlp, target_parameters=p[1]) .. GENERATED FROM PYTHON SOURCE LINES 149-154 The symbolic representation of the parametric sensitivities (jacobian and hessian) is obtained via the :meth:`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``. .. GENERATED FROM PYTHON SOURCE LINES 154-157 .. code-block:: Python J, H = nlp_.parametric_sensitivity(expr=Z, second_order=True) .. GENERATED FROM PYTHON SOURCE LINES 158-161 Once the symbols are available, we can evaluate them numerically after a solution is computed. We do so for different values of :math:`p_2`, and store the results in some lists for later plotting. .. GENERATED FROM PYTHON SOURCE LINES 161-169 .. code-block:: Python 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()) .. GENERATED FROM PYTHON SOURCE LINES 170-176 Visualization of the sensitivities ---------------------------------- Before jumping into the visualization of the sensitivities themselves, let's compute the function :math:`Z` for a range of realizations of :math:`p_2` (with :math:`p_1=0.2`). This will also serve as a baseline to understand whether the sensitivities follow the actual plot. .. GENERATED FROM PYTHON SOURCE LINES 176-181 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 182-187 We can now plot the function :math:`Z` for the different values of :math:`p_2`. After that, we can plot the parametric sensitivities. To do so, we plot the second order Taylor expansion that approximates the function :math:`Z` around the optimal solution for different values of :math:`p_2`. This approximation obviously makes use of the gradient and hessian information of the function at that point. .. GENERATED FROM PYTHON SOURCE LINES 187-207 .. code-block:: Python 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() .. image-sg:: /auto_examples/sensitivity/images/sphx_glr_simple_sensitivity_002.png :alt: simple sensitivity :srcset: /auto_examples/sensitivity/images/sphx_glr_simple_sensitivity_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.055 seconds) .. _sphx_glr_download_auto_examples_sensitivity_simple_sensitivity.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: simple_sensitivity.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: simple_sensitivity.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: simple_sensitivity.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_