Simple markov system get Nan for the first step

Hi

I’m working on an optimal control problem using the acados_template in Python, and I’ve run into an issue that I’m struggling to resolve. The problem is formulated to control an Markov transition system model with constraints on the control input.

Here are some details about my setup:

• The control input is constrained within a certain range.
• I’m using an AcadosOcpSolver with the sqp solver option.
• Even with no constrains the system will get Nan at first step

Unfortunately, when I run my solver, I’m encountering an error with acados returned status 1. I’ve tried to troubleshoot the issue by checking constraints, providing an initial guess, and experimenting with different solver options, but the issue persists.

Here is a snippet of my code where I set up and solve the OCP:

def calculate_motion_mode(temp):
#decompose the rate_matrix
eigenvalues, eigenvectors = np.linalg.eig(get_rate_matrix(temp))

def get_rate_matrix(temp):
E = DM([0, 0.4, 1, 0.2])
B = DM([[inf, 1.5, 1.1, inf],
[1.5, inf, 10, 0.01],
[1.1, 10, inf, 1],
[inf, 0.01, 1, inf]])
rate_matrix = exp(-(B-E)/temp).T

# eliminate the inf values
for i in range(4):
for j in range(4):
if B[i, j] == inf:
rate_matrix[i, j] = 0
for j in range(4):
rate_matrix[j, j] = - sum1(rate_matrix[:, j])
return rate_matrix

def get_equilibrium(temp):
rate_matrix = get_rate_matrix(temp)
# rate 2 transition
trans_matrix = rate_matrix/np.max(np.abs(rate_matrix)) + np.eye(4)
# calculate equilibrium
eigenvalues, eigenvectors = np.linalg.eig(trans_matrix)
eigenvectors = np.transpose(eigenvectors)
targetvector = eigenvectors[np.argmax(eigenvalues)]
targetvector = targetvector / np.sum(targetvector)
return targetvector

model_name = 'asm'
# set up states & controls
x1 = SX.sym('x1', 4)
x = vertcat(x1)

T = SX.sym('T')

# xdot
x1_dot = SX.sym('x1_dot', 4)
xdot = vertcat(x1_dot)

# dynamics
f_expl = vertcat(get_rate_matrix(T)@x1)
f_impl = xdot - f_expl

model.f_impl_expr = f_impl
model.f_expl_expr = f_expl
model.x = x
model.xdot = xdot
model.u = T
model.name = model_name

# set x0 and x_goal
x0 = get_equilibrium(1)
model.x0 = x0
x_goal = get_equilibrium(2)
model.x_goal = x_goal
# motion_mode = calculate_motion_mode(2)
# b = model.x - x_goal
# A = motion_mode[1]
# c = cs.solve(cs.mtimes(A.T, A), cs.mtimes(A.T, b))

# model.con_h_expr_e = vertcat(c[0], c[1])
# model.cost_expr_ext_cost_e = c[2]**2

# store meta information
model.u_labels = ['$T$']
model.t_label = '$t$ [s]'

return model


the solver definition

def main():
Tf = 0.0001

# set model
model = export_asm_model()
ocp.model = model

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

# set number of shooting intervals
ocp.dims.N = N

# set prediction horizon
ocp.solver_options.tf = Tf

# set cost
ocp.constraints.lbu = np.array([1])
ocp.constraints.ubu = np.array([2])
ocp.constraints.idxbu = np.array([0])
ocp.constraints.x0 = model.x0

# set options
ocp.solver_options.qp_solver = 'PARTIAL_CONDENSING_HPIPM' # FULL_CONDENSING_QPOASES
# PARTIAL_CONDENSING_HPIPM, FULL_CONDENSING_QPOASES, FULL_CONDENSING_HPIPM,
# PARTIAL_CONDENSING_QPDUNES, PARTIAL_CONDENSING_OSQP, FULL_CONDENSING_DAQP
ocp.solver_options.hessian_approx = 'GAUSS_NEWTON' # 'GAUSS_NEWTON', 'EXACT'
ocp.solver_options.integrator_type = 'IRK'
# ocp.solver_options.print_level = 1
ocp.solver_options.nlp_solver_type = 'SQP' # SQP_RTI, SQP
ocp.solver_options.globalization = 'MERIT_BACKTRACKING' # turns on globalization

simX = np.zeros((N+1, nx))
simU = np.zeros((N, nu))

status = ocp_solver.solve()


error message

iter    res_stat        res_eq          res_ineq        res_comp        qp_stat qp_iter alpha
0       nan     nan     2.000000e+00    0.000000e+00    0       0       0.000000e+00

Traceback (most recent call last):
File "asm_env/asm_ocp.py", line 68, in <module>
main()
File "asm_env/asm_ocp.py", line 54, in main


Here are my specific questions:

1. What does status code 1 usually indicate in acados?
2. Are there any common pitfalls or typical issues that could lead to this error?
3. Could someone suggest strategies for further debugging or point out potential issues in the code snippet provided?

I would greatly appreciate any advice or insights from those of you who have experience with acados or similar optimal control problems. If additional information is needed, I am happy to provide it.

Hi,

I think the problems are:

1. You did not define any cost. You should define a cost function in order to use acados to solve an optimal control problem.

2. The Nan is probably caused by the fact that you have \exp(-(B-E)/T) in your dynamics. When T is initialized with zeros, this will cause NaNs. So you can try to initialize your T as some non-zeros values.

BTW, since your system dynamics is in the form of \dot x= g(u) x, it would be easier to define \hat u = g(u) and solve for \dot x = \hat u x in acados, then solve u = g^{-1}(\hat u). I believe this is a better way because you get simpler dynamics and you can avoid getting NaNs.