Optimization based solvers

[ErikQQY]

Part of #336

Some TODOs:

• Multi-points BVP are encountering some DI issues, but two-point BVP are fine, need to dig deeper to see what's going wrong.
• Need some patches in SciMLBase.jl and DiffEqBase.jl for better interface.
• More docstrings and documentations

[github-actions]

Benchmark Results

Click to check benchmark results

	master	ab8ea4a...	master / ab8ea4a...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	0.591 ± 0.019 s	0.596 ± 0.014 s	0.992 ± 0.04
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	13.8 ± 0.42 ms	13.8 ± 0.67 ms	1 ± 0.057
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	3.15 ± 0.15 ms	3.13 ± 0.21 ms	1.01 ± 0.082
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	9.37 ± 1.2 ms	9.3 ± 1.3 ms	1.01 ± 0.19
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	1.64 ± 0.28 ms	1.63 ± 0.29 ms	1.01 ± 0.25
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.84 ± 0.65 ms	1.87 ± 0.64 ms	0.986 ± 0.48
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.1 ± 0.97 ms	3.14 ± 0.98 ms	0.986 ± 0.44
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0672 ± 0.0069 s	0.0681 ± 0.0071 s	0.987 ± 0.14
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0802 ± 0.017 s	0.0814 ± 0.0065 s	0.986 ± 0.22
Simple Pendulum/IIP/Shooting(Tsit5())	0.236 ± 0.075 ms	0.245 ± 0.077 ms	0.964 ± 0.43
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.759 ± 0.0046 s	0.759 ± 0.0042 s	1 ± 0.0083
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	16.6 ± 5.7 ms	16.7 ± 5.6 ms	0.999 ± 0.48
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	3.65 ± 0.18 ms	3.62 ± 0.18 ms	1.01 ± 0.072
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	10.9 ± 0.73 ms	10.8 ± 0.76 ms	1.01 ± 0.098
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	1.87 ± 0.14 ms	1.89 ± 0.18 ms	0.99 ± 0.12
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.64 ± 2.9 ms	3.58 ± 2.9 ms	1.02 ± 1.2
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	6.02 ± 5.1 ms	6.07 ± 5.2 ms	0.992 ± 1.2
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.119 ± 0.0043 s	0.12 ± 0.008 s	0.993 ± 0.075
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.145 ± 0.0077 s	0.146 ± 0.006 s	0.995 ± 0.067
Simple Pendulum/OOP/Shooting(Tsit5())	0.641 ± 0.044 ms	0.636 ± 0.047 ms	1.01 ± 0.1
time_to_load	5.4 ± 0.0061 s	5.42 ± 0.061 s	0.996 ± 0.011

### Benchmark Plots A plot of the benchmark results has been uploaded as an artifact to the workflow run for this PR. Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ErikQQY]

The current status of this PR:

Utilization of equality constraints in Optimization.jl

This is just by assigning equal lcons and ucons as zeros so that solvers like Ipopt.jl can be used for the solution of the nonlinear system.✅(Quite similar to what we are doing with NonlinearSolve.jl)

Utilization of inequality constraints in Optimization.jl

This part can cover lots of things:

1. Specify the inequality at a time point. For example, $u_1(\pi/4)<1.0, u_2(\pi/2)>2.0$:✅

function bc!(du, u, p, t)
	res[1] = u(pi/2)[1]
	res[1] = u(pi/4)[2]
end
prob = BVProblem(f!, bc!, u0, tspan, ucons = [1.0, Inf], lcons = [-Inf, 2.0])

2. Specify the whole interval solution to be smaller or bigger than a number: For example, $u_1(t)<1.0$, $u_2(t)>1.0$:✅

function bc!(du, u, p, t)
	res[1] = maxsol(u, (0.0, pi/2))[1]
	res[1] = minsol(u, (0.0, pi/2))[2]
end
prob = BVProblem(f!, bc!, u0, tspan, ucons = [1.0, Inf], lcons = [-Inf, 1.0])

3. Specify the whole interval solution to be smaller or bigger than a continuous function: For example, $u_1(t)<sin(t)$, $u_2(t)>cos(t)$.❌(This functionality is way more complicated than I thought, and involve the reconstruction of the current boundary condition evaluation routines, so might take more time to land this functionality.)

@ChrisRackauckas please review and spot anything I am missing here.