How to use the casadi.SX function in casadi

def shiftfirstknot_T(basis, t_shift, inverse=False):
    # Create transformation matrix that shifts the first (degree+1) knots over
    # t_shift. With inverse = True, the inverse transformation is also
    # computed.
    knots, deg = basis.knots, basis.degree
    N = len(basis)
    if isinstance(t_shift, SX):
        typ, sym = SX, True
    elif isinstance(t_shift, MX):
        typ, sym = MX, True
    else:
        typ, sym = np, False
    _T = typ.eye(deg+1)
    for k in range(deg+1):
        _t = typ.zeros((deg+1+k+1, deg+1+k))
        for j in range(deg+1+k+1):
            if j >= deg+1:
                _t[j, j-1] = 1.
            elif j <= k:
                _t[j, j] = 1.
            else:
                _t[j, j-1] = (knots[j+deg-k]-t_shift)/(knots[j+deg-k]-knots[j])
                _t[j, j] = (t_shift-knots[j])/(knots[j+deg-k]-knots[j])
        _T = mtimes(_t, _T) if sym else _t.dot(_T)

# Model Parameters
    mu1_M = 0.4
    mu2_M = 0.5
    K = 0.05
    K_I = 0.02
    y1 = 0.2
    y2 = 0.15
    Sf = 2.0
    p = 0.5

    # Declare variables (use scalar graph)
    #u =     SX.sym("u")      # parameters

    # Differential states
    xd   = ca.SX.sym("xd", ndstate)  # vector of states [c1,c2,I]

    # Algebraic states
    xa = ca.SX.sym("xa", nastate)  # vector of algebraic states [mu1,mu2]

    # Initial conditions
    xDi = x0

    xAi = [mu1_M * (Sf - x0[0] / y1 - x0[1] / y2) / (K + (Sf - x0[0] / y1 - x0[1] / y2)),
           mu2_M * (Sf - x0[0] / y1 - x0[1] / y2) * K / ((K + (Sf - x0[0] / y1 - x0[1] / y2)) * (K_I + x0[2]))]

    ode = ca.vertcat(xa[0] * xd[0] - xd[0] * u[0],
                     xa[1] * xd[1] - xd[1] * u[0],
                     -p * xd[0] * xd[2] + u[1] - xd[2] * u[0])

    alg = ca.vertcat(mu1_M * (Sf - xd[0] / y1 - xd[1] / y2) / (K + (Sf - xd[0] / y1 - xd[1] / y2)) - xa[0],
                     mu2_M * (Sf - xd[0] / y1 - xd[1] / y2) * K / ((K + (Sf - xd[0] / y1 - xd[1] / y2)) * (K_I + xd[2])) - xa[1])

from acados import ocp_nlp_function, ocp_nlp_ls_cost, ocp_nlp, ocp_nlp_solver
from models import chen_model

N = 10
ode_fun, nx, nu = chen_model(implicit=True)
nlp = ocp_nlp({'N': N, 'nx': nx, 'nu': nu, 'ng': N*[1] + [0]})

# ODE Model
step = 0.1
nlp.set_model(ode_fun, step)

# Cost function
Q = diag([1.0, 1.0])
R = 1e-2
x = SX.sym('x', nx)
u = SX.sym('u', nu)
u_N = SX.sym('u', 0)
f = ocp_nlp_function(Function('ls_cost', [x, u], [vertcat(x, u)]))
f_N = ocp_nlp_function(Function('ls_cost_N', [x, u_N], [x]))
ls_cost = ocp_nlp_ls_cost(N, N*[f]+[f_N])
ls_cost.ls_cost_matrix = N*[block_diag(Q, R)] + [Q]
nlp.set_cost(ls_cost)

# Constraints
g = ocp_nlp_function(Function('path_constraint', [x, u], [u]))
g_N = ocp_nlp_function(Function('path_constraintN', [x, u], [SX([])]))
nlp.set_path_constraints(N*[g] + [g_N])
for i in range(N):
    nlp.lg[i] = -0.5
    nlp.ug[i] = +0.5

solver = ocp_nlp_solver('sqp', nlp, {'sim_solver': 'lifted-irk', 'integrator_steps': 2, 'qp_solver':'hpipm', 'sensitivity_method': 'gauss-newton'})

def create_integrator_cvodes():
    # Time 
    t = C.SX("t")

    # Differential states
    s = C.SX("s"); v = C.SX("v");
    y = [s,v]
    
    # Control
    u = C.SX("u")

    # Parameters
    k = C.SX("k")
    b = C.SX("b")
    m = 1.0   # mass

    # Differential equation - explicit form
    sdot = v
    vdot = (u - k*s - b*v)/m

""" Adjust hard boundries """
        for t in range(Nt):
            for i in range(Ny):
                # Lower boundry of diagonal
                self.__varlb['covariance', t, i + i*Ny] = 0
            self.__varlb['eps', t] = 0
            if xub is not None:
                self.__varub['mean', t] = xub
            if xlb is not None:
                self.__varlb['mean', t] = xlb


        """ Input covariance matrix """
        covar_x_s = ca.MX.sym('covar_x', Ny, Ny)
        covar_x_sx = ca.SX.sym('cov_x', Ny, Ny)
        K_sx = ca.SX.sym('K', Nu, Ny)

        covar_u_func = ca.Function('cov_u', [covar_x_sx, K_sx],
                                   [K_sx @ covar_x_sx @ K_sx.T])

        cov_xu_func = ca.Function('cov_xu', [covar_x_sx, K_sx],
                                  [covar_x_sx @ K_sx.T])

        covar_s = ca.MX(Nx, Nx)
        cov_xu = cov_xu_func(covar_x_s, K_s.reshape((Nu, Ny)))
        covar_s[:Ny, :Ny] = covar_x_s
        covar_s[Ny:, Ny:] = covar_u_func(covar_x_s, K_s.reshape((Nu, Ny)))
        covar_s[Ny:, :Ny] = cov_xu.T
        covar_s[:Ny, Ny:] = cov_xu
        covar_func = ca.Function('covar', [covar_x_s], [covar_s])

n, D = ca.MX.size(X)

    ell = hyper[:D]
    sf2 = hyper[D]**2
    lik = hyper[D + 1]**2
    m   = hyper[D + 2]
    K   = ca.MX.zeros(n, n)

    for i in range(n):
        for j in range(n):
            K[i, j] = covSEard2(X[i, :] - m, X[j, :] - m, ell, sf2)

    K = K + lik * ca.MX.eye(n)
    K = (K + K.T) * 0.5   # Make sure matrix is symmentric

    A = ca.SX.sym('A', ca.MX.size(K))
    #L = ca.chol(A)      # Should check for PD !!!
    cholesky = ca.Function('Cholesky', [A], [ca.chol(A)])
    L = cholesky(K).T

    logK = 2 * ca.sum1(ca.MX.log(ca.fabs(ca.diag(L))))

    invLy = ca.solve(L, Y - m)
    alpha = ca.solve(L.T, invLy)
    NLL = 0.5 * ca.mtimes(Y.T - m, alpha) + 0.5 * logK + n * 0.5 * ca.log(2 * ca.pi)   # - log_priors
    return NLL

def shiftfirstknot_T(basis, t_shift, inverse=False):
    # Create transformation matrix that shifts the first (degree+1) knots over
    # t_shift. With inverse = True, the inverse transformation is also
    # computed.
    knots, deg = basis.knots, basis.degree
    N = len(basis)
    if isinstance(t_shift, SX):
        typ, sym = SX, True
    elif isinstance(t_shift, MX):
        typ, sym = MX, True
    else:
        typ, sym = np, False
    _T = typ.eye(deg+1)
    for k in range(deg+1):
        _t = typ.zeros((deg+1+k+1, deg+1+k))
        for j in range(deg+1+k+1):
            if j >= deg+1:
                _t[j, j-1] = 1.
            elif j <= k:
                _t[j, j] = 1.
            else:
                _t[j, j-1] = (knots[j+deg-k]-t_shift)/(knots[j+deg-k]-knots[j])
                _t[j, j] = (t_shift-knots[j])/(knots[j+deg-k]-knots[j])

def sx_sym(name, dim1 = 1, dim2 = 1):

    return ca.SX.sym(name, dim1, dim2)

def shiftoverknot(self):
    T, knots = self.getBasis().get_shiftoverknot_tf()
    coeffs = self.getCoefficient().getData()
    if isinstance(coeffs, (SX, MX)):
        coeffs2 = mtimes(T, coeffs)
    else:
        coeffs2 = T.dot(coeffs)
    basis2 = spl.BSplineBasis(knots, self.getBasis().getDegree())
    return spl.Function(basis2, coeffs2)

z += ocp.p
   z += ocp.t
   g += list(F.eval([z])[0])
   for j in range(n_x):
       g.append(dxs[i][j] - (xs[i+1][j] - xs[i][j])/(tf/N))

cost = xs[N][0]

# Equality constraint residual
g_res = casadi.SXFunction([xx],[g])

# Objective function
obj = casadi.SXFunction([xx], [[cost]])

# Hessian
sigma = casadi.SX('sigma')

lam = []
Lag = sigma*cost
for i in range(len(g)):
   lam.append(casadi.SX('lambda_' + str(i)))
   Lag = Lag + g[i]*lam[i]

Lag_fcn = casadi.SXFunction([xx, lam, [sigma]],[[Lag]])
H = Lag_fcn.hess()

#H_fcn = casadi.SXFunction([xx, lam, [sigma]],[H])
H_fcn = Lag_fcn.hessian(0,0)

# Create Solver
solver = casadi.IpoptSolver(obj,g_res,H_fcn)