odow
commented
9 months ago I can't run your code because I don't know what PCY1 is.
As a general comment, JuMP v0.18.6 was released over four years ago. In the last few years we have made big improvements to our nonlinear code, so that the expression and eval stuff is not needed. See https://jump.dev/JuMP.jl/stable/manual/nonlinear/ and tutorials like https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/rocket_control/. I should write an explicit MPC tutorial (x-ref https://github.com/jump-dev/JuMP.jl/issues/2348#issuecomment-1839210372).
You can ask for support and suggested code improvements on http://jump.dev/forum. We try to keep the GitHub issues for bug reports.
Here's how I'd write your code. I can't test it, but it should point you in the right direction:
using JuMP, Ipopt
# TODO start
PCY1 = 1.0
PDY1 = 1.0
PDY2 = 1.0
PEY1 = 1.0
PEY2 = 1.0
PEY3 = 1.0
PKY1 = 1.0
PKY2 = 1.0
PHY1 = 1.0
PHY2 = 1.0
PVY1 = 1.0
PVY2 = 1.0
FZ0 = 1.0
EP_SMALL = 1e-3
# TODO end
KZX = 289.5
KZYR = 427.0
KZYF = 277.2
FzF0 = 2670.0
FzR0 = 2670.0
PC1 = PCY1
PD1 = PDY1 - PDY2
PD2 = PDY2/FZ0
PE1 = PEY1 - PEY2
PE2 = PEY2/FZ0
PE3 = PEY3
PK1 = PKY1*FZ0
PK2 = 1/(PKY2*FZ0)
PH1 = PHY1 - PHY2
PH2 = PHY2/FZ0
PV1 = PVY1 - PVY2
PV2 = PVY2/FZ0
g = 9.81
T = 0.002
la = 1.5521
lb = 2.715 - la #1.1629
h_cg = 0.596
m = 1.269e+03
Izz = 1.620e+03
FzF0 = 2670.0
FzR0 = 2670.0
L = 3.6
TF = 4
num_states = 8
num_controls = 2
state_name = [:xpos, :ypos, :v, :r, :psi, :sa, :ux, :ax]
control_name = [:sr, :jx]
XL = [-10, -3, -3, -pi, -pi, -pi/4, 0.01, -3]
XU = [200, 4, 3, pi, pi, pi/4, 10, 3]
X0 = [0, 0, 0, 0, 0.001, 0, 2, 0]
CL = [-pi/4 * 5, -2.6]
CU = [pi/4 * 5, 2.6]
n = 25
Δt = TF / n
function ThreeDOFBicycle_expr(x::Vector, u::Vector)
xpos, ypos, v, r, psi, sa, ux, ax = x
sr, jx = u
FYF = (
PD2 * (FzF0 - (ax - v*r)*KZX)^2 +
PD1 * (FzF0 - (ax - v*r)*KZX)
) * sin(
PC1 * atan(
(((PK1*sin(2*atan(PK2*(FzF0 - (ax - v*r)*KZX))))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)) + ((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX))^2 + EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v + la*r)/(ux+EP_SMALL)) - sa) + PH2*(FzF0 - (ax - v*r)*KZX) + PH1)) - ((PE2*(FzF0 - (ax - v*r)*KZX) + PE1)*(1 - PE3)*(((atan((v + la*r)/(ux+EP_SMALL)) - sa) + PH2*(FzF0 - (ax - v*r)*KZX) + PH1))/((((atan((v + la*r)/(ux+EP_SMALL)) - sa) + PH2*(FzF0 - (ax - v*r)*KZX) + PH1)^2 + EP_SMALL^2)^(0.5)))*((((PK1*sin(2*atan(PK2*(FzF0 - (ax - v*r)*KZX))))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)) + ((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX))^2 + EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v + la*r)/(ux+EP_SMALL)) - sa) + PH2*(FzF0 - (ax - v*r)*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzF0 - (ax - v*r)*KZX))))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)) + ((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX)))/(((PD2*PC1*(FzF0 - (ax - v*r)*KZX)^2 + PD1*PC1*(FzF0 - (ax - v*r)*KZX))^2 + EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v + la*r)/(ux+EP_SMALL)) - sa) + PH2*(FzF0 - (ax - v*r)*KZX) + PH1)))))) + (PV2*(FzF0 - (ax - v*r)*KZX)^2 + PV1*(FzF0 - (ax - v*r)*KZX))
FYR = (
PD2 * (FzR0 + (ax - v*r)*KZX)^2 +
PD1 * (FzR0 + (ax - v*r)*KZX)
) * sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzR0 + (ax - v*r)*KZX))))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)) + ((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX))^2+EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v - lb*r)/(ux+EP_SMALL))) + PH2*(FzR0 + (ax - v*r)*KZX) + PH1)) - ((PE2*(FzR0 + (ax - v*r)*KZX) + PE1)*(1 - PE3*(((atan((v - lb*r)/(ux+EP_SMALL))) + PH2*(FzR0 + (ax - v*r)*KZX) + PH1))/((((atan((v - lb*r)/(ux+EP_SMALL))) + PH2*(FzR0 + (ax - v*r)*KZX) + PH1)^2 + EP_SMALL^2)^(0.5))))*((((PK1*sin(2*atan(PK2*(FzR0 + (ax - v*r)*KZX))))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)) + ((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX))^2+EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v - lb*r)/(ux+EP_SMALL))) + PH2*(FzR0 + (ax - v*r)*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzR0 + (ax - v*r)*KZX))))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)) + ((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX)))/(((PD2*PC1*(FzR0 + (ax - v*r)*KZX)^2 + PD1*PC1*(FzR0 + (ax - v*r)*KZX))^2+EP_SMALL^2)^(0.5))*0.001)+EP_SMALL))*((atan((v - lb*r)/(ux+EP_SMALL))) + PH2*(FzR0 + (ax - v*r)*KZX) + PH1)))))) + (PV2*(FzR0 + (ax - v*r)*KZX)^2 + PV1*(FzR0 + (ax - v*r)*KZX))
return [
ux * cos(psi) - (v + la * r) * sin(psi),
ux * sin(psi) + (v + la * r) * cos(psi),
(FYF + FYR) / m - r * ux,
(la * FYF - lb * FYR) / Izz,
r,
sr,
ax,
jx,
]
end
simp_veh = Model(Ipopt.Optimizer)
set_silent(simp_veh)
@variables(simp_veh, begin
XL[i] <= x[j in 1:n, i in 1:num_states] <= XU[i]
CL[i] <= u[j in 1:n, i in 1:num_controls] <= CU[i]
end)
fix.(x[1, :], X0; force = true)
dx = [ThreeDOFBicycle_expr(x[j, :], u[j, :]) for j in 1:n]
@constraint(
simp_veh,
[st in 1:num_states, j in 1:n-1],
x[j+1,st] == x[j,st] + Δt * dx[j+1][st],
)
@objective(
simp_veh,
Min,
5*(x[n,1] - 100)^2 + sum(Δt * (x[i,2] - 3.6)^2 + 10 * Δt * u[i,1]^2 + Δt * x[i,3]^2 for i in 1:n)^2,
)
println("Solving...")
optimize!(simp_veh)
println("Min Obj: ", solve_time(simp_veh))
for i in 1:1:10
optimize!(simp_veh)
println("Min Obj: ", solve_time(simp_veh))
endNote that the final for-loop isn't going to be interesting as a comparison to v0.18.6, because Ipopt will now detect that the problem structure hasn't changed and so it won't re-build the problem.
If you can provide a reproducible example, there is still plenty of room for improvement, because you use a lot of non-constant global variables.
johnysyumich
Hi developers,
We are currently "new" users to JuMP but we do have experience in using JuMP since we were users of NLOptControl: https://github.com/JuliaMPC/NLOptControl.jl. This package relies very heavily on JuMP and stopped upgrading about three years ago. We are thinking of upgrading it to let it compile with the current version of JuMP. The last version of NLOptControl uses JuMP v0.18.6 which is the version right before the change from MPB to MOI. However, we met some performance issues when we were upgrading it, and here is one of the dummy examples we created to benchmark with the old package.
Here is the code we run: example.zip
Basically, it is controlling a modeled vehicle to change to another lane.
On my computer (I7-12700K), the algorithm can solve the problem with around 130ms. However, with the old implementation (NLOptControl & JuMP v0.18.6) we can achieve around 45ms per optimization on the same computer. With our current modifications, we cannot achieve the same performance as we could before. It would be very helpful if you would provide any insights about the difference between performance and how we can modify the current dummy code to work as fast as the previous version.
Thanks for your help in advance!