Multiplication of matrix expression and variables leads to stack overflow and matmul error

[GNCGenie]

I have an optimal orbital transfer problem:

h = 2
n = 6
m = 3

model = Model(NLopt.Optimizer)
set_optimizer_attribute(model, "algorithm", :AUGLAG)
local_optimizer = NLopt.Opt(:LD_LBFGS, n*(h+1) + m*h)
local_optimizer.xtol_rel = 1e-4
set_optimizer_attribute(model, "local_optimizer", local_optimizer)

@variables model begin
X[1:6, 1:h+1]
U[1:3, 1:h]
27000.0 <= T[1:h]
end

Trying to multiply an expression matrix

B
6×3 Matrix{Any}:
  ((2 X[1,1]²) / (X[1,1] * sqrt(3.986004415e14 / X[1,1]))) * (X[2,1] * sin(X[6,1]))                 …  0
  (1.0 / (X[1,1] * sqrt(3.986004415e14 / X[1,1]))) * ((X[1,1] * (-X[2,1]² + 1)) * sin(X[6,1]))         0
 0                                                                                                      (((X[1,1] * (-X[2,1]² + 1)) / (1.0 + (X[2,1] * cos(X[6,1])))) * cos(X[6,1] + X[4,1])) / (X[1,1] * sqrt(3.986004415e14 / X[1,1]))
  (-(X[1,1] * (-X[2,1]² + 1)) * cos(X[6,1])) / (X[2,1] * (X[1,1] * sqrt(3.986004415e14 / X[1,1])))      -((((X[1,1] * (-X[2,1]² + 1)) / (1.0 + (X[2,1] * cos(X[6,1])))) * sin(X[6,1] + X[4,1])) * cos(X[3,1])) / ((X[1,1] * sqrt(3.986004415e14 / X[1,1])) * sin(X[3,1]))
 0                                                                                                      (((X[1,1] * (-X[2,1]² + 1)) / (1.0 + (X[2,1] * cos(X[6,1])))) * sin(X[6,1] + X[4,1])) / ((X[1,1] * sqrt(3.986004415e14 / X[1,1])) * sin(X[3,1]))
  ((X[1,1] * (-X[2,1]² + 1)) * cos(X[6,1])) / (X[2,1] * (X[1,1] * sqrt(3.986004415e14 / X[1,1])))   …  0

with

U[:, i] = JuMP.VariableRef[U[1,1], U[2,1], U[3,1]]

by executing:

B1 = B(X[:,1])
B1*U[:,1]

Results in the following stack overflow error:

operate(::typeof(*), ::Matrix{Any}, ::Vector{JuMP.VariableRef})@LinearAlgebra.jl:403
mul(::Matrix{Any}, ::Vector{JuMP.VariableRef})@shortcuts.jl:58
*(::Matrix{Any}, ::Vector{JuMP.VariableRef})@dispatch.jl:360
 [ Repeated 100 times ]
operate(::typeof(*), ::Matrix{Any}, ::Vector{JuMP.VariableRef})@LinearAlgebra.jl:403
mul(::Matrix{Any}, ::Vector{JuMP.VariableRef})@shortcuts.jl:58
*(::Matrix{Any}, ::Vector{JuMP.VariableRef})@dispatch.jl:360

Whereas doing:

B1*U[:,:]

Gives the correct result, a 6x2 matrix.

Help in understanding and resolving the problem is appreciated!

[jd-foster]

It's hard to reproduce your issue since your code is also not reproducible: can you either give the full expression for B or a function that produces it? At one point where you write B1 = B(X[:,1]) it also seems that B is a function? It is also probably not necessary to have an Optimizer set for the model as this is working at the algebraic level.

[GNCGenie]

Definitely! The code to reproduce the issue:

    using JuMP
    using NLopt

    const R_EARTH     = 6.378136300e6              # [m] GGM05s Value
    const GM_EARTH    = 3.986004415e14          # [m^3/s^2] GGM05s Value
    const J2_EARTH    = 0.0010826358191967      # [] GGM05s value
    const OMEGA_EARTH = 7.292115146706979e-5    
    ωₑ = OMEGA_EARTH;
    μ = GM_EARTH;
    Rₑ = R_EARTH;
    J₂ = J2_EARTH;

function dT(x...)
	(a,e,i,ω,Ω,θ) = x
	p = a*(1-e^2)
	r = p/(1+e*cos(θ))
	h = a*√(μ/a)
	ωₛ = √(μ/a^3)

	A = replace(
		[.0,.0,.0,.0,(-3/2)*(Rₑ^2/(a*(1-e^2))^2)*J₂*ωₛ*cos(i),h/r^2]
		, NaN=>0, Inf=>0, -Inf=>0)
	
	return A
end

function dU(x...)
	(a,e,i,ω,Ω,θ) = x
	p = a*(1-e^2)
	r = p/(1+e*cos(θ))
	h = a*√(μ/a)
	ωₛ = √(μ/a^3)
	
	B = replace(
		[(2*a^2/h)*(e*sin(θ)) (2*a^2/h)*(p/r) 0
		(1/h)*(p*sin(θ)) (1/h)*((p+r)*cos(θ) + r*e) 0
		0 0 (r*cos(θ+ω))/h
		-p*cos(θ)/(e*h) (p+r)*sin(θ)/(e*h) -(r*sin(θ+ω)*cos(i))/(h*sin(i))
		0 0 (r*sin(θ+ω))/(h*sin(i))
		(p*cos(θ))/(e*h) -((p+r)*sin(θ))/(e*h) 0]
		, NaN=>0, Inf=>0, -Inf=>0)
	
	return B
end

begin
	h = 2
	n = 6
	m = 3
	
	model = Model(NLopt.Optimizer)
	set_optimizer_attribute(model, "algorithm", :AUGLAG)
	local_optimizer = NLopt.Opt(:LD_LBFGS, n*(h+1) + m*h)
	local_optimizer.xtol_rel = 1e-4
	set_optimizer_attribute(model, "local_optimizer", local_optimizer)

	@variables model begin
		X[1:6, 1:h+1]
		U[1:3, 1:h]
		27000.0 <= T[1:h]
	end
	let i=1
                B = dU(X[:,i]...)
		B*U[:,i]
	end
end

This is required to set up the dynamics constraints for the control problem, U is the effort input and B is the jacobian of X with effort. Later on when the formulation is done, we minimise the norm of U to minimise fuel use of transfer.

[jd-foster]

I can run your code, but I can't reproduce any errors on the latest JuMP version. I recommend upgrading to the latest package versions. On Julia Version 1.10.2, these are

(jl_6E3tDY) pkg> st
Status `../jl_6E3tDY/Project.toml`
  [4076af6c] JuMP v1.20.0
  [76087f3c] NLopt v1.0.2

[GNCGenie]

Yes! That resolved the issue.