Updating parameters after the OCP solver has been created has no effect

Hi :wave:

I use acados with the python interface and my issue refers to updating parameters after the solver is created:

I have noticed that the parameter values set via ocp.parameter_values before the OCP solver is created effects the solution although I updated the parameters for all stages after the solver was created and before solve() was called.

The phenomenon can be reproduced using the test_parametric_nonlinear_constaint_h.py script

Setting initial values to p_0

ocp.parameter_values = p_0
# create solver
p_set = p_0 #np.zeros(n_param)
for i in range(N):
    ocp_solver.set(i, "p", p_set)
# solve()

yield the following result graph:

Setting initial values to np.zeros(n_param and later updating with correct values p_0

ocp.parameter_values = np.zeros(n_param)
# create solver
p_set = p_0 
for i in range(N):
    ocp_solver.set(i, "p", p_set)
# solve()

yield slightly different results:

Any help is much appreciated, thanks!

Attached the full code to test the issue (run in folder ./acados/examples/acados_python/tests):

import sys
sys.path.insert(0, '../pendulum_on_cart/common')

from acados_template import AcadosOcp, AcadosOcpSolver
from pendulum_model import export_pendulum_ode_model
import numpy as np
import scipy.linalg
from utils import plot_pendulum
from casadi import SX, vertcat

def main():
    ocp = AcadosOcp()

    # set model
    model = export_pendulum_ode_model()
    ocp.model = model
    x = model.x
    u = model.u

    Tf = 1.0
    nx = x.rows()
    nu = u.rows()
    N = 20

    # set dimensions
    ocp.solver_options.N_horizon = N

    # set cost
    Q = 2*np.diag([1e3, 1e3, 1e-2, 1e-2])
    R = 2*np.diag([1e-2])
    cost_W = scipy.linalg.block_diag(Q, R)

    #
    n_param = 42
    p = SX.sym('p', n_param)
    constraint_quotient = p[0]
    y_param = p[1:nx+nu+1]
    ocp.model.p = p

    # define cost with parametric reference
    ocp.cost.cost_type = 'EXTERNAL'
    ocp.cost.cost_type_e = 'EXTERNAL'

    residual = y_param - vertcat(x, u)
    ocp.model.cost_expr_ext_cost = residual.T @ cost_W @ residual
    res_e = y_param[0:nx] - x
    ocp.model.cost_expr_ext_cost_e = res_e.T @ Q @ res_e

    # set constraints

    Fmax = 80

    ocp.constraints.lh = np.array([-Fmax])
    ocp.constraints.uh = np.array([+Fmax])
    ocp.model.con_h_expr = model.u / constraint_quotient

    ocp.constraints.lh_0 = np.array([-Fmax])
    ocp.constraints.uh_0 = np.array([+Fmax])
    ocp.model.con_h_expr_0 = model.u / constraint_quotient

    p_0 = np.zeros(n_param)
    p_0[0] = 1.0
    p_0[1] = 1.0

    ### CHANGE THIS TO TEST THE BUG ###
    # ocp.parameter_values = p_0
    ocp.parameter_values = np.zeros(n_param)

    ocp.constraints.x0 = np.array([0.0, np.pi, 0.0, 0.0])

    ocp.solver_options.qp_solver = 'PARTIAL_CONDENSING_HPIPM' # FULL_CONDENSING_QPOASES
    ocp.solver_options.hessian_approx = 'EXACT' # GAUSS_NEWTON, EXACT
    ocp.solver_options.regularize_method = 'CONVEXIFY' # GAUSS_NEWTON, EXACT
    ocp.solver_options.integrator_type = 'ERK'
    ocp.solver_options.print_level = 0
    ocp.solver_options.nlp_solver_type = 'SQP' # SQP_RTI, SQP

    # set prediction horizon
    ocp.solver_options.tf = Tf

    # Cython
    if 0:
        AcadosOcpSolver.generate(ocp, json_file='acados_ocp.json')
        AcadosOcpSolver.build(ocp.code_export_directory, with_cython=True)
        ocp_solver = AcadosOcpSolver.create_cython_solver('acados_ocp.json')
    else:
        ocp_solver = AcadosOcpSolver(ocp, json_file = 'acados_ocp.json')
    

    print("--- update params ---")
    p_set = p_0 #np.zeros(n_param)
    for i in range(N):
        ocp_solver.set(i, "p", p_set)

    solX = np.zeros((N+1, nx))
    solU = np.zeros((N, nu))

    print("--- solve ocp ---")
    status = ocp_solver.solve()

    if status != 0:
        raise Exception(f'acados returned status {status}.')

    # get solution
    for i in range(N):
        solX[i,:] = ocp_solver.get(i, "x")
        solU[i,:] = ocp_solver.get(i, "u")
    solX[N,:] = ocp_solver.get(N, "x")

    # ocp_solver.print_statistics()

    if np.any(solU > Fmax) or np.any(solU < -Fmax):
        raise Exception(f"control bounds should be respected by solution, got u: {solU}")

    if not np.allclose(np.max(solU), Fmax):
        raise Exception(f"control should go to bounds, got u: {solU}")
    if not np.allclose(np.min(solU), -Fmax):
        raise Exception(f"control should go to bounds, got u: {solU}")

    plot_pendulum(np.linspace(0, Tf, N+1), Fmax, solU, solX, latexify=False)


if __name__ == "__main__":
    main()

I think I found my issue:
There are N_horizon + 1 parameter stages.

1 Like