Undeclared identifier DISCRETE in generated C code

Hi there,

I am trying to implement an ocp to control the CSTR system. I would like to use discrete dynamics, but the generated C code does not compile and returns the following error:

acados_sim_solver_cstr.c:86:23: error: use of undeclared identifier 'DISCRETE'
    plan.sim_solver = DISCRETE;
                      ^
1 error generated.
make: *** [acados_sim_solver_cstr.o] Error 1

You can find below my whole implementation. I thought that specifying disc_dyn_expr in the model and integrator_type = “DISCRETE” in the ocp instance would suffice, especially since all the example files (including getting_started/example_ocp_dynamics_formulations.py) run normally.

Thank you a lot for your input!

Tudor

from acados_template import *
from casadi import *
from scipy.linalg import solve_discrete_are

xr1 = np.array(
    [
        2.1402105301746182e00,
        1.0903043613077321e00,
        1.1419108442079495e02,
        1.1290659291045561e02,
    ]
)
ur1 = np.array([14.19, -1113.50])
xr2 = np.array([2.7151681, 1.02349152, 105.50585058, 100.8920758])
ur2 = np.array([13.66640639, -3999.58908628])
xr3 = np.array([2.97496989, 0.95384459, 101.14965441, 95.19386292])
ur3 = np.array([12.980148, -5162.95653026])


def export_cstr_model(dt: float, rk4_nodes: int = 6):
    c_A = SX.sym("c_A")
    c_B = SX.sym("c_B")
    theta = SX.sym("theta")
    theta_K = SX.sym("theta_K")
    u_1 = SX.sym("u_1")
    u_2 = SX.sym("u_2")
    x = vertcat(c_A, c_B, theta, theta_K)
    u = vertcat(u_1, u_2)
    c_A_dot = SX.sym("c_A_dot")
    c_B_dot = SX.sym("c_B_dot")
    theta_dot = SX.sym("theta_dot")
    theta_K_dot = SX.sym("theta_K_dot")
    xdot = vertcat(c_A_dot, c_B_dot, theta_dot, theta_K_dot)

    k_10 = 1.287e12  # [h^-1]
    k_20 = 1.287e12  # [h^-1]
    k_30 = 9.043e9  # [h^-1]
    E_1 = -9758.3  # [1]
    E_2 = -9758.3  # [1]
    E_3 = -8560.0  # [1]
    H_1 = 4.2  # [kJ.mol^-1]
    H_2 = -11.0  # [kJ.mol^-1]
    H_3 = -41.85  # [kJ.mol^-1]
    rho = 0.9342  # [kg.L^-1]
    C_p = 3.01  # [kJ/(kg.K)]
    k_w = 4032.0  # [kJ/(h.m^2.K)]
    A_R = 0.215  # [m^2]
    V_R = 10.0  # [l]
    m_K = 5.0  # [kg]
    C_PK = 2.0  # [kJ/(kg.K)]
    c_A0 = 5.1  # [mol/l]
    theta_0 = 104.9  # [°C]

    k_1 = k_10 * exp(E_1 / (x[2] + 273.15))
    k_2 = k_20 * exp(E_2 / (x[2] + 273.15))
    k_3 = k_30 * exp(E_3 / (x[2] + 273.15))

    f_cont = Function(
        "f_cont",
        [x, u],
        [
            vertcat(
                u[0] * (c_A0 - x[0]) - k_1 * x[0] - k_3 * x[0] * x[0],
                -u[0] * x[1] + k_1 * x[0] - k_2 * x[1],
                u[0] * (theta_0 - x[2])
                + (x[3] - x[2]) * k_w * A_R / (rho * C_p * V_R)
                - (k_1 * x[0] * H_1 + k_2 * x[1] * H_2 + k_3 * x[0] * x[0] * H_3)
                / (rho * C_p),
                (u[1] + k_w * A_R * (x[2] - x[3])) / (m_K * C_PK),
            )
        ],
    )
    new_x = x
    for j in range(rk4_nodes):
        k1 = f_cont(new_x, u)
        k2 = f_cont(new_x + dt / 2 * k1, u)
        k3 = f_cont(new_x + dt / 2 * k2, u)
        k4 = f_cont(new_x + dt * k3, u)
        new_x += dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

    model = AcadosModel()
    model.f_impl_expr = xdot - f_cont(x, u)
    model.f_expl_expr = f_cont(x, u)
    model.disc_dyn_expr = new_x
    model.x = x
    model.xdot = xdot
    model.u = u
    model.name = "cstr"

    return model


def export_cstr_ocp(
    dt: float = 20 / 3600,
    N: int = 100,
    x0: np.ndarray = np.array([1.0, 0.5, 100.0, 100.0]),
    x_ref: np.ndarray = xr1,
    u_ref: np.ndarray = ur1,
):
    ocp = AcadosOcp()
    ocp.model = export_cstr_model(dt=dt)
    ocp.dims.N = N
    ocp.solver_options.tf = dt * N

    nx = ocp.model.x.size()[0]
    nu = ocp.model.u.size()[0]
    x = ocp.model.x
    u = ocp.model.u

    Q = np.diag([0.2, 1.0, 0.5, 0.2])
    R = np.diag([0.5, 5.0e-7])
    C_x = np.vstack((np.eye(nx), -np.eye(nx)))
    d_x = np.array([10.0, 10.0, 150.0, 150.0, 0.0, 0.0, -98.0, -92.0])
    q_x = 2 * nx
    C_u = np.vstack((np.eye(nu), -np.eye(nu)))
    d_u = np.array([35.0, 0.0, -3.0, 9000.0])
    q_u = 2 * nu
    epsilon = 1.0e-3

    # cost function
    jac_f_x = Function("jac_f_x", [x, u], [jacobian(ocp.model.disc_dyn_expr, x)])
    jac_f_u = Function("jac_f_u", [x, u], [jacobian(ocp.model.disc_dyn_expr, u)])
    A = np.array(jac_f_x(x_ref, u_ref))
    B = np.array(jac_f_u(x_ref, u_ref))
    ocp.cost.cost_type = "LINEAR_LS"
    ocp.cost.cost_type_e = "LINEAR_LS"
    ocp.cost.W = np.block([[Q, np.zeros((nx, nu))], [np.zeros((nu, nx)), R]])
    P = np.array(solve_discrete_are(A, B, Q, R))
    ocp.cost.W_e = P
    ocp.cost.Vx = np.vstack((np.eye(nx), np.zeros((nu, nx))))
    ocp.cost.Vx_e = np.eye(nx)
    ocp.cost.Vu = np.vstack((np.zeros((nx, nu)), np.eye(nu)))
    ocp.cost.yref = np.append(x_ref, u_ref)
    ocp.cost.yref_e = x_ref

    # constraints
    ocp.constraints.x0 = x0
    ocp.constraints.lbx = d_x[:nx]
    ocp.constraints.ubx = d_x[nx:]
    ocp.constraints.idxbx = np.arange(nx)
    ocp.constraints.lbu = d_u[:nu]
    ocp.constraints.ubu = d_u[nu:]
    ocp.constraints.idxbu = np.arange(nu)

    # solver options
    # ocp.solver_options.qp_solver = "FULL_CONDENSING_QPOASES"
    ocp.solver_options.qp_solver = "PARTIAL_CONDENSING_HPIPM"
    ocp.solver_options.hessian_approx = "GAUSS_NEWTON"
    # ocp.solver_options.integrator_type = "ERK"
    ocp.solver_options.integrator_type = "DISCRETE"
    ocp.solver_options.nlp_solver_type = "SQP_RTI"

    return ocp


def run_closed_loop_simulation(
    dt: float = 20 / 3600,
    N: int = 100,
    x0: np.ndarray = np.array([1.0, 0.5, 100.0, 100.0]),
    xr: np.ndarray = xr1,
    ur: np.ndarray = ur1,
):
    ocp = export_cstr_ocp(dt=dt, N=N, x0=x0, x_ref=xr, u_ref=ur)
    acados_ocp_solver = AcadosOcpSolver(
        ocp, json_file="acados_ocp_" + ocp.model.name + ".json"
    )
    acados_sim_solver = AcadosSimSolver(
        ocp, json_file="acados_sim_" + ocp.model.name + ".json"
    )

    Nsim = 350
    nx = ocp.model.x.size()[0]
    nu = ocp.model.u.size()[0]
    x_sim = np.zeros((Nsim + 1, nx))
    u_sim = np.zeros((Nsim, nu))
    x_sim[0, :] = x0

    xcurrent = x0

    n_convergence = Nsim + 1
    for i in range(Nsim):
        # solve ocp
        acados_ocp_solver.set(0, "lbx", xcurrent)
        acados_ocp_solver.set(0, "ubx", xcurrent)
        status = acados_ocp_solver.solve()
        acados_ocp_solver.print_statistics()
        if status != 0:
            raise Exception(
                "acados ocp solver returned status {}. Exiting.".format(status)
            )

        u_sim[i, :] = acados_ocp_solver.get(0, "u")

        # simulate system
        acados_sim_solver.set("x", xcurrent)
        acados_sim_solver.set("u", u_sim[i, :])
        status = acados_sim_solver.solve()
        if status != 0:
            raise Exception(
                "acados sim solver returned status {}. Exiting.".format(status)
            )

        # update state
        xcurrent = acados_sim_solver.get("x")
        x_sim[i + 1, :] = xcurrent

        # check if there is convergence in relative norm
        if np.linalg.norm(x_sim[i + 1, :] - xr) / np.linalg.norm(xr) < 1e-3:
            n_convergence = i + 1
            break

    x_sim = x_sim[: n_convergence + 1, :]
    u_sim = u_sim[:n_convergence, :]

    return x_sim, u_sim, n_convergence


if __name__ == "__main__":
    run_closed_loop_simulation()

Hi Tudor,

the DISCRETE dynamics module is not supported within the AcadosSimSolver class.
I think that your example fails in this command:

    acados_sim_solver = AcadosSimSolver(
        ocp, json_file="acados_sim_" + ocp.model.name + ".json"
    )

You can just use the discrete dynamics function directly in your closed loop simulation instead of an acados integrator.

Best,
Jonathan

That’s what I did and it worked, thank you very much.

Best,
Tudor