Hello,
I have been trying to improve my solve time for nonlinear optimal control problem , it is a tool for solving minimum lap time for a race car GitHub - ceserik/SLapSim.jl: Scalable laptime simulator · GitHub. I am using JuMP + IPOPT + MUMPS, I have got licence for HSL and pardiso solver, but then I realized most of the time is not spent solving the linear system of equations by linear solver, but computing the hessian of my function. So the different linear solvers had no impact on time of solving.
This is a snippet from my code, where I enforce gauss-lobatto quadrature for states along time.
#create constraints using butcher tableau of lobattoIIIA
segment_start_idx = 1
# evaluate the car ODE function and adds state specific constraints for first point
f(X[segment_start_idx, :], U[segment_start_idx, :], s_all[segment_start_idx], model)
for segment = 1:number_of_segments
F = [f(X[segment_start_idx+col-1, :], U[segment_start_idx+col-1, :], s_all[segment_start_idx+col-1], nothing) for col in 1:stages]
for stage = 2:stages
fxSum = zeros(NonlinearExpr, nStates)
for col = 1:stages
fxSum += tableau.a[stage, col] * F[col]
end
# enforce loabbto quadrature
@constraint(model, X[segment_start_idx+stage-1, :] .== X[segment_start_idx, :] + h_all[segment] * fxSum)
# evaluate the car ODE function and adds state specific constraints for rest of the points
f(X[segment_start_idx+stage-1, :], U[segment_start_idx+stage-1, :], s_all[segment_start_idx+stage-1], model)
end
segment_start_idx = segment_start_idx + stages - 1
end
This snippet is of tire function used by car, it constraints the tire force according to complex tire model using trigionometry functions
function compute(inTorque::carVar, optiModel::Union{JuMP.Model,Nothing}=nothing)
slipAngle.value = -atan(velocity.value[2], velocity.value[1])
α = slipAngle.value * 1
D = forces.value[3] * frictionCoefficient.value
#Pacejka tire model
forces.value[2] = D * sin(C * atan(B * α - E * (B * α - atan(B * α))))
forces.value[1] = inTorque/radius.value + brakingForce.value + forces.value[3] * rollingResistance.value
end
function tireConstraints(model=nothing)
μFz = frictionCoefficient.value * forces.value[3]
fx = forces.value[1] / μFz
fy = forces.value[2] / μFz
@constraint(model, fx^2 + fy^2 <= 1.1^2)
end
Time statistics of solver
Number of Iterations....: 318
(scaled) (unscaled)
Objective...............: 4.3202620034505081e+01 8.6405240069010162e+01
Dual infeasibility......: 4.4689843001849362e-09 8.9379686003698724e-09
Constraint violation....: 2.6980963047834181e-09 1.1728004665201297e-06
Variable bound violation: 9.2138580376488477e-09 9.2138580376488477e-09
Complementarity.........: 1.0103392529261279e-09 2.0206785058522557e-09
Overall NLP error.......: 4.4689843001849362e-09 1.1728004665201297e-06
Number of objective function evaluations = 391
Number of objective gradient evaluations = 319
Number of equality constraint evaluations = 391
Number of inequality constraint evaluations = 391
Number of equality constraint Jacobian evaluations = 319
Number of inequality constraint Jacobian evaluations = 319
Number of Lagrangian Hessian evaluations = 318
Total seconds in IPOPT (w/o function evaluations) = 10.715
Total seconds in NLP function evaluations = 166.565
Timing Statistics:
OverallAlgorithm....................: 235.910 (sys: 1.131 wall: 177.280)
PrintProblemStatistics.............: 0.000 (sys: 0.000 wall: 0.000)
InitializeIterates.................: 0.287 (sys: 0.004 wall: 0.291)
UpdateHessian......................: 121.763 (sys: 0.059 wall: 121.749)
OutputIteration....................: 0.015 (sys: 0.002 wall: 0.012)
UpdateBarrierParameter.............: 0.138 (sys: 0.004 wall: 0.084)
ComputeSearchDirection.............: 38.595 (sys: 0.981 wall: 9.960)
ComputeAcceptableTrialPoint........: 67.673 (sys: 0.056 wall: 37.725)
AcceptTrialPoint...................: 0.030 (sys: 0.000 wall: 0.030)
CheckConvergence...................: 7.400 (sys: 0.025 wall: 7.421)
PDSystemSolverTotal.................: 38.502 (sys: 0.979 wall: 9.942)
PDSystemSolverSolveOnce............: 36.372 (sys: 0.966 wall: 9.413)
ComputeResiduals...................: 2.013 (sys: 0.010 wall: 0.495)
StdAugSystemSolverMultiSolve.......: 36.175 (sys: 0.962 wall: 9.387)
LinearSystemScaling................: 0.000 (sys: 0.000 wall: 0.000)
LinearSystemSymbolicFactorization..: 0.031 (sys: 0.002 wall: 0.033)
LinearSystemFactorization..........: 0.000 (sys: 0.000 wall: 0.000)
LinearSystemBackSolve..............: 4.002 (sys: 0.006 wall: 1.032)
LinearSystemStructureConverter.....: 0.000 (sys: 0.000 wall: 0.000)
LinearSystemStructureConverterInit: 0.000 (sys: 0.000 wall: 0.000)
QualityFunctionSearch...............: 0.000 (sys: 0.000 wall: 0.000)
TryCorrector........................: 0.000 (sys: 0.000 wall: 0.000)
Task1...............................: 0.000 (sys: 0.000 wall: 0.000)
Task2...............................: 0.000 (sys: 0.000 wall: 0.000)
Task3...............................: 0.000 (sys: 0.000 wall: 0.000)
Task4...............................: 0.000 (sys: 0.000 wall: 0.000)
Task5...............................: 0.000 (sys: 0.000 wall: 0.000)
Task6...............................: 0.000 (sys: 0.000 wall: 0.000)
Function Evaluations................: 196.333 (sys: 0.134 wall: 166.565)
Objective function.................: 0.006 (sys: 0.000 wall: 0.006)
Objective function gradient........: 0.018 (sys: 0.001 wall: 0.018)
Equality constraints...............: 67.240 (sys: 0.054 wall: 37.462)
Inequality constraints.............: 0.137 (sys: 0.000 wall: 0.138)
Equality constraint Jacobian.......: 0.039 (sys: 0.000 wall: 0.039)
Inequality constraint Jacobian.....: 7.175 (sys: 0.020 wall: 7.193)
Lagrangian Hessian.................: 121.718 (sys: 0.058 wall: 121.708)
EXIT: Optimal Solution Found.
I want to ask if it is normal that 90% of ipopt run time is spent in function evaluation( Hessian and equality constraints in the bottom) and very little time solving the linear system by HSL/MUMPS/PARDISO. Does this hint I formulated the model improperly?
If it is useful for somebody hessian and jacobian look like this
