Hello, I am making a tool to find fastest possible lap time for racing car, I am trying to use same function for initial value problem solving and also for optimisation in julia. The function describes behavior of a car. I have a rotation matrix to rotate a wheel, steering angle is therefore union of float and optimization variable. I am getting this error when I added steering of wheels to the optimisation.
[1] operate(::typeof(*), A::Matrix{NonlinearExpression}, B::Vector{Union{Float64, NonlinearExpression, VariableRef, AffExpr}})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/PcwmY/src/implementations/LinearAlgebra.jl:427
[2] mul(::Matrix{NonlinearExpression}, ::Vector{Union{Float64, NonlinearExpression, VariableRef, AffExpr}})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/PcwmY/src/shortcuts.jl:58
[3] *(a::Matrix{NonlinearExpression}, b::Vector{Union{Float64, NonlinearExpression, VariableRef, AffExpr}})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/PcwmY/src/dispatch.jl:369--- the above 3 lines are repeated 26660 more times ---
and this is the minimal working example that causes the same error.
Hi @Ceserik, you’re using the legacy @NL API. This does not support linear algebra. Use the new API, which is just removing @NL, and changing NonlinearExpression to JuMP.NonlinearExpr.
Hello, thanks for your answer. In my actual code I am using new API, I created minimal example using legacy API by mistake. And I forgot to use actual rotation matrix. I am sorry for not describing my issue better. I have written new minimal example, that is closer to what I am actually running. The rotation works fine for Float64 but not for JuMP expressions.
using JuMP
model = JuMP.Model()
v = Vector{Union{Float64, JuMP.VariableRef, JuMP.AffExpr, JuMP.NonlinearExpression}}(undef, 2)
w = Vector{Union{Float64, JuMP.VariableRef, JuMP.AffExpr, JuMP.NonlinearExpression}}(undef, 2)
@variable(model,x[1:2])
@variable(model, angle)
v[1] = x[1]
v[2] = x[2]
w[1] = 1.0
w[2] = 2.0
function rotate(angle)
out = [
cos(angle) -sin(angle);
-sin(angle) cos(angle)
]
return out
end
b = rotate(pi/2) * w #works fine
c = rotate(angle) * v #error
By the way this function responsible for creating and solving the problem, I cannot include it all here because the car model is spread around many files
function findOptimalTrajectory(track::Track,car::Car,model::JuMP.Model)
N = length(track.sampleDistances)
#fill track parameters which are constant along track, should be automatized
track.rho = fill(track.rho[1],N)
track.μ = fill(track.μ[1],N)
track.widthL = fill(track.widthL[1],N)
track.widthR = fill(track.widthR[1],N)
nStates = Int(car.carParameters.nStates.value)
#create inputs for car model #create instance of track parameters
#determine sizes of inputs and states
nControls = Int(round(car.carParameters.nControls.value))
@variable(model,U[1:N-1,1:nControls])
#stateInit = [fill(5,N),zeros(N),zeros(N),zeros(N),zeros(N),range(0,4,N)]'
@variable(model,X[1:N,1:6])
for i = 1:N
set_start_value(X[i,1], 5.0) # vx = 5
set_start_value(X[i,2], 0.0) # vy = 0
set_start_value(X[i,3], pi/2) # psi = 0
set_start_value(X[i,4], 0.0) # dpsi = 0
set_start_value(X[i,5], 0.0) # n = 0
set_start_value(X[i,6], 4*(i-1)/(N-1)) # t = linspace(0,4,N)
#println(4*(i-1)/(N-1))
end
s = track.sampleDistances#1:length(track.curvature)
#t_final = X[end,6]
for k = 1:N-1 # loop over control intervals
x_next = X[k,:] .+ (s[k+1]-s[k]) .* carODE_path(car,track,k, U[k,:], X[k,:],model);
@constraint(model,X[k+1,:] .== x_next)
end
#initial states
# [vx vy psi dpsi n t]
#initial = [4; 0; 0; 0; 0; 0]
#set initial states
#@constraint(model,X[1,:] .== initial)
@constraint(model,X[1:end,1] .>= 0) #vx
@constraint(model,X[1,1] .== 5) # intial vx
@constraint(model,X[1,2] .== 0) # intial vx
@constraint(model,X[1,6] .>= 0) # final time
@constraint(model,diff(X[:,6]) .>=0) #time goes forward
@objective(model,Min,X[end,6])
optimize!(model)
# Plotting the results
fig = Figure(layout = GridLayout(6, 1))
# Create subplots for each variable
labels = ["vx", "vy", "psi", "dpsi", "n", "t"]
for (i, label) in enumerate(labels)
ax = Axis(fig[i, 1], ylabel = label)
lines!(ax, value.(X[:,i]), label=label)
end
display(GLMakie.Screen(), fig) # This creates a new window
return X, U
end
