Comparison of CasADi’s and csnlp’s sensitivity computations

In this example, we reproduce the sensitivity computations from this CasADi example, both with casadi’s native capabilities as well as the one from csnlp.

This example assumes you are already familiar with the basics of computing sensitivities via csnlp.wrappers.NlpSensitivity. If this is not the case, please refer to the example A simple example of sensitivity analysis for a more introductory example.

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

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.

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()

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.

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)

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.

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)

-------------------------------------------
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

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.

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()

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

With 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.

nlp.init_solver({"print_time": False, "ipopt": {"print_level": 0, "sb": "yes"}})
sol_ = nlp.solve(pars={"p": 0})

Then, we can compute the sensitivities of the solution w.r.t. the parameter p with the csnlp.wrappers.NlpSensitivity class. By passing a solution object to the csnlp.wrappers.NlpSensitivity.parametric_sensitivity method, computations are sped up because they are numerical.

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()

Comparison

We can now compare the results of the two methods.

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",
)

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]]

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.

print(
    "\nTimings", "CasADi time:", f" {t1:.6f} s", "csnlp time:", f" {t3:.6f} s", sep="\n"
)

Timings
CasADi time:
 0.009239 s
csnlp time:
 39.774833 s