.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/sensitivity/comparison_sensitivities.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_comparison_sensitivities.py: Comparison of CasADi's and :mod:`csnlp`'s sensitivity computations ================================================================== In this example, we reproduce the sensitivity computations from `this CasADi example `_, both with :mod:`casadi`'s native capabilities as well as the one from :mod:`csnlp`. This example assumes you are already familiar with the basics of computing sensitivities via :class:`csnlp.wrappers.NlpSensitivity`. If this is not the case, please refer to the example :ref:`simple_sensitivity_example` for a more introductory example. .. GENERATED FROM PYTHON SOURCE LINES 13-23 .. code-block:: Python import time import casadi as cs import numpy as np from numpy.polynomial import Polynomial from csnlp import Nlp, wrappers from csnlp.core.data import array2cs, cs2array .. GENERATED FROM PYTHON SOURCE LINES 24-28 With CasADi's native functionality ---------------------------------- We start by defining the collocation coefficients, the dynamics and cost functions, and some constants of the optimal control problem. .. GENERATED FROM PYTHON SOURCE LINES 28-72 .. code-block:: Python def get_coeffs(d: int) -> tuple[cs.DM, cs.DM, cs.DM]: """Gets the coefficients of the collocation equation, of the continuity equation, and lastly of the quadrature function.""" tau_root = np.append(0, cs.collocation_points(d, "legendre")) C = cs.DM.zeros(d + 1, d + 1) D = cs.DM.zeros(d + 1) B = cs.DM.zeros(d + 1) for j in range(d + 1): p = Polynomial([1]) for r in range(d + 1): if r != j: p *= Polynomial([-tau_root[r], 1]) / (tau_root[j] - tau_root[r]) D[j] = p(1.0) pder = p.deriv() for r in range(d + 1): C[j, r] = pder(tau_root[r]) pint = p.integ() B[j] = pint(1.0) return C, D, B def get_dynamics_and_cost() -> cs.Function: """Returns the dynamics and objective function of the problem.""" x1 = cs.SX.sym("x1") x2 = cs.SX.sym("x2") x = cs.vertcat(x1, x2) u = cs.SX.sym("u") xdot = cs.vertcat((1 - x2**2) * x1 - x2 + u, x1) L = x1**2 + x2**2 + u**2 return cs.Function("f", [x, u], [xdot, L], ["x", "u"], ["xdot", "L"]) # time and control horizon T = 10.0 N = 6 # number of control intervals h = T / N # coefficients and dynamics and cost functions d = 3 C, D, B = get_coeffs(d) F = get_dynamics_and_cost() .. GENERATED FROM PYTHON SOURCE LINES 73-76 Then, we build the problem itself. This time we use collocation points to enforce the dynamics and the continuity of the states. We also add a perturbation parameter to the initial state. .. GENERATED FROM PYTHON SOURCE LINES 76-123 .. code-block:: Python nlp = Nlp[cs.SX]() # create variables nx, nu = 2, 1 X, _, _ = nlp.variable("x", (nx, N + 1), lb=[[-0.25], [-np.inf]]) XC, _, _ = nlp.variable("x_colloc", (nx, N * d), lb=[[-0.25], [-np.inf]]) U, _, _ = nlp.variable("u", (nu, N), lb=-1, ub=0.85) # create initial state and perturbation parameter nlp.constraint("init", X[:, 0], "==", [0, 1]) p = nlp.parameter("p", (nx, 1)) # formulate NLP with collocation J = 0.0 X_ = cs2array(X).T XC_ = cs2array(XC).reshape((nx, N, d)).transpose((1, 0, 2)) U_ = cs2array(U).T for k, (x, x_next, xc, u) in enumerate(zip(X_[:-1], X_[1:], XC_, U_)): # convert back to symbolic x = array2cs(x) x_next = array2cs(x_next) xc = array2cs(xc) u = array2cs(u) # perturb first state with p if k == 0: x += p # propagate collocation points and compute associated cost f_k, q_k = F(xc, u) # create collocation constraints xp = x @ C[0, 1:] + xc @ C[1:, 1:] nlp.constraint(f"colloc_{k}", h * f_k, "==", xp) # Add contribution to quadrature function J += h * (q_k @ B[1:]) # create end state and constraint x_k_end = D[0] * x + xc @ D[1:] nlp.constraint(f"colloc_end_{k}", x_next, "==", x_k_end) # set objective nlp.minimize(J) .. GENERATED FROM PYTHON SOURCE LINES 124-130 Now that the problem has been created, we can use the high-level approach to sensitivity, native to CasADi. To use it, it is mandatory to use the ``sqpmethod`` solver. First, we convert the problem to a ``sqpmethod``-based problem, and solve it numerically. .. GENERATED FROM PYTHON SOURCE LINES 130-142 .. code-block:: Python prob = {"f": nlp.f, "x": nlp.x, "g": nlp.g, "p": nlp.p} opts = { "qpsol": "qrqp", "qpsol_options": {"print_iter": False, "error_on_fail": False}, "print_time": False, } solver = cs.nlpsol("solver", "sqpmethod", prob, opts) kwargs = {"lbx": nlp.lbx.data, "ubx": nlp.ubx.data, "lbg": 0, "ubg": 0, "p": 0} sol = solver(x0=0, **kwargs) .. rst-class:: sphx-glr-script-out .. code-block:: none ------------------------------------------- This is casadi::QRQP Number of variables: 56 Number of constraints: 50 Number of nonzeros in H: 78 Number of nonzeros in A: 260 Number of nonzeros in KKT: 662 Number of nonzeros in QR(V): 693 Number of nonzeros in QR(R): 1051 ------------------------------------------- This is casadi::Sqpmethod. Using exact Hessian Number of variables: 56 Number of constraints: 50 Number of nonzeros in constraint Jacobian: 260 Number of nonzeros in Lagrangian Hessian: 78 iter objective inf_pr inf_du ||d|| lg(rg) ls info 0 0.000000e+00 1.00e+00 0.00e+00 0.00e+00 - 0 - 1 3.832471e+00 2.00e-01 3.05e+00 8.00e-01 - 2 - WARNING(sqpmethod): Indefinite Hessian detected 2 4.887941e+00 4.00e-02 7.67e-01 2.00e-01 - 1 - 3 4.783708e+00 3.49e-05 1.18e-03 1.96e-02 - 1 - 4 4.783682e+00 1.97e-10 2.29e-09 3.37e-05 - 1 - MESSAGE(sqpmethod): Convergence achieved after 4 iterations .. GENERATED FROM PYTHON SOURCE LINES 143-146 We then spawn, with the native factory method, a new solver that yields the jacobian and hessian of the objective w.r.t. NLP parameters. Since we have already compute a solution, we can use it to warm-start the new solver. .. GENERATED FROM PYTHON SOURCE LINES 146-163 .. code-block:: Python t0 = time.perf_counter() # create the jacobian and hessian solver jac_and_hess_solver = solver.factory("s", solver.name_in(), ["jac:f:p", "hess:f:p:p"]) # call it with the previous solution to warm-start it jac_and_hess_sol = jac_and_hess_solver( x0=sol["x"], lam_x0=sol["lam_x"], lam_g0=sol["lam_g"], **kwargs ) t1 = time.perf_counter() - t0 # retrieve the arrays jac_f_p = jac_and_hess_sol["jac_f_p"].full().squeeze() hess_f_p_p = jac_and_hess_sol["hess_f_p_p"].full() .. rst-class:: sphx-glr-script-out .. code-block:: none iter objective inf_pr inf_du ||d|| lg(rg) ls info 0 4.783682e+00 1.97e-10 2.29e-09 0.00e+00 - 0 - MESSAGE(sqpmethod): Convergence achieved after 0 iterations .. GENERATED FROM PYTHON SOURCE LINES 164-171 With :mod:`csnlp` ---------------------------------- We'll now turn our attention on how to compute the same sensitivities with our package. Since it does not rely on the ``sqpmethod`` solver, we can use any, especially IPOPT. First, let us initialize IPOPT and compute the same solution. .. GENERATED FROM PYTHON SOURCE LINES 171-175 .. code-block:: Python nlp.init_solver({"print_time": False, "ipopt": {"print_level": 0, "sb": "yes"}}) sol_ = nlp.solve(pars={"p": 0}) .. GENERATED FROM PYTHON SOURCE LINES 176-180 Then, we can compute the sensitivities of the solution w.r.t. the parameter ``p`` with the :class:`csnlp.wrappers.NlpSensitivity` class. By passing a solution object to the :meth:`csnlp.wrappers.NlpSensitivity.parametric_sensitivity` method, computations are sped up because they are numerical. .. GENERATED FROM PYTHON SOURCE LINES 180-192 .. code-block:: Python nlp_ = wrappers.NlpSensitivity(nlp) t2 = time.perf_counter() dfdp, d2fdp2 = nlp_.parametric_sensitivity(nlp.f, solution=sol_, second_order=True) t3 = time.perf_counter() - t2 dfdp = dfdp.full().squeeze() d2fdp2 = d2fdp2.full() .. GENERATED FROM PYTHON SOURCE LINES 193-196 Comparison ---------- We can now compare the results of the two methods. .. GENERATED FROM PYTHON SOURCE LINES 196-206 .. code-block:: Python print( "\nComparison", "Abs. differences in Jacobian of f w.r.t. p:", np.abs(jac_f_p.T - dfdp), "Abs. differences in Hessian of f w.r.t. p:", np.abs(hess_f_p_p - d2fdp2), sep="\n", ) .. rst-class:: sphx-glr-script-out .. code-block:: none Comparison Abs. differences in Jacobian of f w.r.t. p: [1.13183241e-10 1.05238485e-10] Abs. differences in Hessian of f w.r.t. p: [[8.93152219e-12 7.05808745e-11] [7.05719927e-11 4.75530726e-12]] .. GENERATED FROM PYTHON SOURCE LINES 207-211 The price we pay for such a general approach is that our method is invariably slower than the native one, especially as the problem size grows with ``N``. Specifically, the computations of the hessian are the slowest. The jacobian computations can sustain larger sizes. .. GENERATED FROM PYTHON SOURCE LINES 211-215 .. code-block:: Python print( "\nTimings", "CasADi time:", f" {t1:.6f} s", "csnlp time:", f" {t3:.6f} s", sep="\n" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Timings CasADi time: 0.009239 s csnlp time: 39.774833 s .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 39.818 seconds) .. _sphx_glr_download_auto_examples_sensitivity_comparison_sensitivities.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: comparison_sensitivities.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: comparison_sensitivities.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: comparison_sensitivities.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_