# The Jacobian for the equality constraints contains an invalid number (+ general help with matrices in JuMP)

Hello all! I am pretty new to Julia, and am trying to implement a matrix form model for the dynamics. This is what I have written so far:

``````using JuMP #, MadNLP
import Ipopt

#Number of mesh point
n = 10

#Simulation parameters
const g = [-1; 0; 0]
const m_wet = 2.0
const m_dry = 1.0
const T_min = 0.3
const T_max = 5
const δ_max = 20
const θ_max = 90
const γ_gs = 20
const ω_max = 60
const J_B = 1e-2 * [1 0 0 ; 0 1 0 ; 0 0 1]
const r_TB = [-1e-2 0 0]

const I_sp = 1                    # s

#Initial and final conditions
r₀ = [4 4 0]
v_f = [-1exp(-1) 0 0]
ω₀ = [0 0 0]
q_f = [1 0 0 0]

#Initial calculations

αₘ = 1 / I_sp / -9.80665

#Defining model
model = Model(Ipopt.Optimizer)

#Defining statez variables
@variable(model, r[1:n, 1:3])
@variable(model, q[1:n, 1:4])
@variable(model, v[1:n, 1:3])
@variable(model, ω[1:n, 1:3])
@variable(model, m_dry ≤ m[1:n] ≤ m_wet)
@variable(model, 1e-4 ≤ t_f)

#Defining control variables
@variable(model, T_B[1:n, 1:3])

@constraint(model, T_min^2 .≤ sum(T_B.^2, dims=2) .≤ T_max^2)

#Dynamic equation expressions
#Matrix Functions
function C_IB_func(q)
C = [(1-2*(q[2]^2+q[3]^2)) 2*(q[1]*q[2] - q[4]*q[3])  2*(q[1]*q[3] + q[4]*q[2]);
2*(q[1]*q[2] + q[4]*q[3]) (1-2*(q[1]^2+q[3]^2)) 2*(q[2]*q[3] - q[4]*q[1]);
2*(q[1]*q[3] - q[4]*q[2]) 2*(q[2]*q[3] + q[4]*q[1]) (1-2*(q[1]^2+q[2]^2))]
return C
end

function Omega(ξ)
A = [0 -ξ[1] -ξ[2] -ξ[3];
ξ[1] 0 ξ[3] -ξ[2];
ξ[2] -ξ[3] 0 ξ[2];
ξ[3] ξ[2] -ξ[1] 0;]

return A
end

@expression(model, Δτ, t_f/(n-1))                       # is this the best way to do that?

@expression(model, δm[j=1:n], -αₘ * sqrt(sum(T_B[j, :].^2)))

function δv_calculations(q, T, g, m)
C = C_IB_func(q)
matrixres = g + C*T./m
return matrixres
end

@expression(model, δv[j=1:n, k=1:3], δv_calculations(q[j,:], T_B[j,:], g, m[j])[k])

function δq_calculations(ω, q)
return 0.5 * Omega(ω) * q
end

@expression(model, δq[j=1:n, k=1:4], δq_calculations(ω[j, :], q[j,:])[k])

function ξcross(ξ)
M = [0 -ξ[3] ξ[2];
ξ[3] 0 -ξ[1];
-ξ[2] ξ[1] 0]
return M
end

function δω_calculations(J_B, r_TB, T, ω)
A = ξcross(r_TB)*T - ξcross(ω)*J_B*ω
B = inv(J_B) * A
return B
end

@expression(model, δω[j=1:n, k=1:3], δω_calculations(J_B, r_TB, T_B[j, :], ω[j, :])[k])

integration_rule = "trapezoidal"

for j=2:n
i = j-1
for k=1:3
@NLconstraint(model, -r[j, k] + r[i, k] + 0.5 * Δτ * (v[j, k] + v[i, k]) == 0)
@NLconstraint(model, -ω[j, k] + ω[i, k] + 0.5 * Δτ * (δω[j, k] + δω[i, k]) == 0)
@NLconstraint(model, -v[j, k] + v[i, k] + 0.5 * Δτ * (δv[j, k] + δv[i, k]) == 0)
end
for k=1:4
@NLconstraint(model, -q[j, k] + q[i, k] + 0.5 * Δτ * (δq[j, k] + δq[i, k]) == 0)
end
@NLconstraint(model, -m[j] + m[i] + 0.5 * Δτ * (δm[j] + δm[i]) == 0)
end

@objective(model, Max, m[n])

optimize!(model)

``````

Even though I don’t get any errors defining the problem, whenever I try to run the last optimization stage, I get the error that:
“The Jacobian for the equality constraints contains an invalid number”

I have read in this same forum that this is normally due to a division by zero, but I couldn’t find any points where I have done this.

I have identified the problem comes from the last non linear constraint (the one related to m). Could it be from including the norm in the previous constraint? Is there a better way to define this?

Also, is the way I have dealt with matrix calculations the best? It is the “cleanest” way I could get it to work in terms of lines of code, but I am not sure if it is the best numerically.

When ipopt starts, everything is set to zero. In the code, in the constraint:

``````    @NLconstraint(model, -m[j] + m[i] + 0.5 * Δτ * (δm[j] + δm[i]) == 0)
``````

if `δm=0`, then `m=0` and the function `δv_calculations` tries to divide by zero, so it fails. Now, `δm=0` if `T_B=0` which it shouldn’t thanks to this constraint:

``````@constraint(model, T_min^2 .≤ sum(T_B.^2, dims=2) .≤ T_max^2)
``````

So when we force `T_B` different from zero in the first iteration by setting

``````@variable(model, T_B[1:n, 1:3],start=1.0) #####Can be something else, doesn't need to be 1.0
``````

then the solution is found:

Also, welcome to the forum

1 Like

Hi @Jack21,

Just to follow-up on what @blob suggests, here’s how I would write your model:

``````using JuMP
import Ipopt
n = 10
const g = [-1; 0; 0]
const m_wet = 2.0
const m_dry = 1.0
const T_min = 0.3
const T_max = 5
const δ_max = 20
const θ_max = 90
const γ_gs = 20
const ω_max = 60
const J_B = 1e-2 * [1 0 0; 0 1 0; 0 0 1]
const r_TB = [-1e-2, 0, 0]
const I_sp = 1
r₀ = [4, 4, 0]
v_f = [-1exp(-1), 0, 0]
ω₀ = [0, 0, 0]
q_f = [1, 0, 0, 0]
αₘ = 1 / I_sp / -9.80665

function C_IB_func(q)
return [
(1-2*(q[2]^2+q[3]^2)) 2*(q[1]*q[2] - q[4]*q[3]) 2*(q[1]*q[3] + q[4]*q[2])
2*(q[1]*q[2] + q[4]*q[3])     (1-2*(q[1]^2+q[3]^2)) 2*(q[2]*q[3] - q[4]*q[1])
2*(q[1]*q[3] - q[4]*q[2]) 2*(q[2]*q[3] + q[4]*q[1])     (1-2*(q[1]^2+q[2]^2))
]
end
function Omega(ξ)
return [
0    -ξ[1] -ξ[2] -ξ[3]
ξ[1]     0  ξ[3] -ξ[2]
ξ[2] -ξ[3]     0  ξ[2]
ξ[3]  ξ[2] -ξ[1]    0
]
end
ξcross(ξ) = [0 -ξ[3] ξ[2]; ξ[3] 0 -ξ[1]; -ξ[2] ξ[1] 0]

model = Model(Ipopt.Optimizer)
@variables(model, begin
r[1:n, 1:3]
q[1:n, 1:4]
v[1:n, 1:3]
ω[1:n, 1:3]
m_dry <= m[1:n] <= m_wet
1e-4 <= t_f
T_B[1:n, 1:3], (start = 1)
end)
@expressions(model, begin
Δτ, t_f / (n-1)
δm[j in 1:n], -αₘ * sqrt(sum(T_B[j, :].^2))
δv[j in 1:n], g .+ C_IB_func(q[j, :]) * T_B[j, :] ./ m[j]
δq[j in 1:n], 0.5 * Omega(ω[j, :]) * q[j, :]
δω[j in 1:n],
inv(J_B) * (ξcross(r_TB)*T_B[j, :] - ξcross(ω[j, :])*J_B*ω[j, :])
end)
@constraints(model, begin
[j in 1:n], T_min^2 <= sum(T_B[j, :].^2) <= T_max^2
[j in 2:n, k in 1:3], r[j, k] == r[j-1, k] + 0.5 * Δτ * (v[j, k] + v[j-1, k])
[j in 2:n, k in 1:3], ω[j, k] == ω[j-1, k] + 0.5 * Δτ * (δω[j][k] + δω[j-1][k])
[j in 2:n, k in 1:3], v[j, k] == v[j-1, k] + 0.5 * Δτ * (δv[j][k] + δv[j-1][k])
[j in 2:n, k in 1:4], q[j, k] == q[j-1, k] + 0.5 * Δτ * (δq[j][k] + δq[j-1][k])
[j in 2:n], m[j] == m[j-1] + 0.5 * Δτ * (δm[j] + δm[j-1])
end)
@objective(model, Max, m[n])
optimize!(model)
``````
1 Like

Dear @blob and @odow,

Thank you so much for your time and help, this is super helpful!