Jacobian Error. Evaluating Jacobian of inequality constraints at user provided starting point. No scaling factors.

[mario35to]

Hello all,

I have an error "WARNING("nlp:nlp_jac_g failed: NaN detected for output jac_g_x, at nonzero index 1 (row 2, col 0).") [.../casadi/core/oracle_function.cpp:377]" when it try to implement MPC optimization of a modelica model .mo file.
For Casadi, I use Pymoca to convert the model in Jupyter Python. I refer the following link to do optimization: https://github.com/pymoca/pymoca/blob/master/test/notebooks/Casadi.ipynb

I couldn't find an answer regarding to this error with my model. Because there is no inconsistent operation in OMEdit model below such as square-root, division by zero etc.
"oscillating masses" model is created in OMEdit.
OMEdit Code:
model oscillating_masses
parameter Real k1=2;
parameter Real k2=2;
parameter Real m1=0.5;
parameter Real m2=1.5;
parameter Real c1=1;
parameter Real c2=1;
input Real u(start = 1);
Real x1 (start = 2);
Real x2 (start = 1.5);
Real x3 (start = 1);
Real x4 (start = 0.5);
initial equation
x1 = 2.0;
x2 = 1.5;
x3 = 1.0;
x4 = 0.5;
equation
der(x1) = x3;
der(x2) = x4;
der(x3) = -((k1+k2)/m1)*x1 + (k2/m1)*x2 -((c1+c2)/m1)*x3 + (c2/m1)*x4 + (1/m1)*u;
der(x4) = (k2/m2)*x1 -(k2/m2)*x2 + (c2/m2)*x3 - (c2/m2)*x4;
end oscillating_masses;

Casadi code:

%matplotlib inline
from pymoca.backends.casadi.api import transfer_model
import numpy as np
from casadi import *
import pylab
from pylab import plot, legend
#load .mo file
model = transfer_model('C:\Users\ekint\AppData\Roaming\Python\Python311\site-packages\pymoca\test\models\', 'oscillating_masses', {})
dae = model.dae_residual_function #connection of dae function and Modelica model
X0 = [model.states[0].start, model.states[1].start, model.states[2].start, model.states[3].start] #initial value

nx = 4 #number of states
nu= 1 #number of input
Q = 20
Q = Qnp.diag(np.ones(nx))
R = 10
R = np.diag(Rnp.ones(nu))
x = SX.sym("x",nx,1)
u = SX.sym("u",nu,1)

Discretization steps

N = 15

Discretization time step

dt = 0.1

Optimization variables

X = MX.sym('X',(N)nx,1)
U = MX.sym('U',Nnu,1)
#state boundaries
lbx = [-4.0, -10.0, -4.0, -10.0]
ubx = [4.0, 10.0, 4.0, 10.0]
#input boundaries
lbu = [-1.5]
ubu = [1.5]
stage_cost = x.T@Q@x+ u.T@R@u
stage_cost_fcn = Function('stage_cost',[x,u],[stage_cost])
lb_g = [] # lower bound for constraint expression g
ub_g = [] # upper bound for constraint expression g
g = [] # constraint expression
f = 0 #initial objective function
lbX = [] # lower bound for X
ubX = [] # Upper bound for X
lbU = [] # lower bound for U
ubU = [] # Upper bound for U
k = 0 #initial value for optimization loop

#initial value for objective function
f = stage_cost_fcn(X[k*nx:(k+1)nx,:],U[knu:(k+1)nu,:])
print(f)
#residual function that collocate DAE using backwards Euler method
res = dae(0, X[knx:(k+1)nx,:], (X[knx:(k+1)nx,:] - X0) / dt, U[0], U[knu:(k+1)*nu,:], MX(), MX())
#appending residual function to constraint
g.append(res)
#appending zeros to lower and upper bound of constraint
lb_g.append(np.zeros((nx,1)))
ub_g.append(np.zeros((nx,1)))
#appending bounds of states and inputs
lbX.append(lbx)
ubX.append(ubx)
lbU.append(lbu)
ubU.append(ubu)

#creating optimization problem with loop
for k in range(1, N):
x_k = X[k*nx:(k+1)nx,:]
u_k = U[knu:(k+1)nu,:]
f += stage_cost_fcn(x_k,u_k) #objective function
print(stage_cost_fcn(x_k,u_k))
res = dae(k * dt, X[knx:(k+1)nx,:], (X[knx:(k+1)*nx,:] - X[(k-1)*nx:(k)nx,:]) / dt, MX(), U[knu:(k+1)*nu,:], MX(), MX())
g.append(res)
lb_g.append(np.zeros((nx,1)))
ub_g.append(np.zeros((nx,1)))
lbX.append(lbx)
ubX.append(ubx)
lbU.append(lbu)
ubU.append(ubu)

lbx = vertcat(*lbX, *lbU)
ubx = vertcat(*ubX, *ubU)
lbg = vertcat(*lb_g)
ubg = vertcat(*ub_g)
nlp = {'x': ca.vertcat(X, U), 'f': f, 'g': ca.vertcat(*g)}
ipopt_options = {'tol': 1e-4} #ipopt(Interior Point OPTimizer) option settings
solver = nlpsol('nlp', 'ipopt', nlp, {'ipopt': ipopt_options})
solution = solver(lbx=lbx, ubx=ubx, lbg=lbg, ubg=ubg)

Full error message: ******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt

---

This is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...: 251
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 75

Error evaluating Jacobian of equality constraints at user provided starting point.
No scaling factors for equality constraints computed!

Number of Iterations....: 0

Number of objective function evaluations = 0
Number of objective gradient evaluations = 0
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 0.001

EXIT: Invalid number in NLP function or derivative detected.
nlp : t_proc (avg) t_wall (avg) n_eval
nlp_grad_f | 0 ( 0) 33.00us ( 33.00us) 1
nlp_jac_g | 1.00ms (500.00us) 798.00us (399.00us) 2
total | 3.00ms ( 3.00ms) 2.48ms ( 2.48ms) 1
CasADi - 2024-03-26 15:56:57 WARNING("nlp:nlp_jac_g failed: NaN detected for output jac_g_x, at nonzero index 1 (row 2, col 0).") [.../casadi/core/oracle_function.cpp:377]
CasADi - 2024-03-26 15:56:57 WARNING("nlp:nlp_jac_g failed: NaN detected for output jac_g_x, at nonzero index 1 (row 2, col 0).") [.../casadi/core/oracle_function.cpp:377]

Any ideas about this error?

Thanks,
Ekin