Optimal control interface

[ErikQQY]

Close: #336

The current implementation concatenates the inequality/equality(the upper part) and the loss function([loss_bc; loss_collocation] as the lower part) as the constraints in OptimizationProblem. The Jacobian is handled by concatenating the Jacobian of inequality constraints(the upper square matrix) and the Jacobian of the loss function ([jac_loss_bc; jac_loss_collocation] as the lower rectangular matrix), resulting in a (N*M + length(bc_resid_prototype) + length(f_prototype)*(N-1)) x (N*M) Jacobian. Then the Jacobian is, upper part: banded matrix; lower part: almost banded matrix.

There are some weird issues with formulating the state variables and control variables, which may prevent the solver from finding an optimal solution. I’ll see what’s going wrong here.

Some TODOs:

• Better ucons and lcons handling. More friendly information on choosing ucons and lcons should be provided to users.
• Better cost function definition. The cost function definition should allow for using the interpolations, e.g. cost(sol, p) = sol(pi)[2]-1.0. Maybe we can also allow integral expressions.

[github-actions]

Benchmark Results

Click to check benchmark results

	master	c02fa42...	master / c02fa42...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	0.32 ± 0.021 s	0.318 ± 0.024 s	1.01 ± 0.1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	9.89 ± 2.1 ms	9.72 ± 2 ms	1.02 ± 0.3
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	2.51 ± 0.14 ms	2.5 ± 0.11 ms	1.01 ± 0.071
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	7.37 ± 1.8 ms	7.3 ± 1.9 ms	1.01 ± 0.36
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	1.33 ± 0.36 ms	1.31 ± 0.37 ms	1.02 ± 0.4
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	2.37 ± 1 ms	1.83 ± 1 ms	1.29 ± 0.91
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.09 ± 1.5 ms	3.03 ± 1.2 ms	1.02 ± 0.63
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0737 ± 0.024 s	0.0736 ± 0.026 s	1 ± 0.48
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.1 ± 0.034 s	0.0986 ± 0.032 s	1.02 ± 0.48
Simple Pendulum/IIP/Shooting(Tsit5())	0.249 ± 0.07 ms	0.25 ± 0.074 ms	0.996 ± 0.41
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.503 ± 0.0096 s	0.492 ± 0.007 s	1.02 ± 0.024
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	12.8 ± 7.4 ms	12.7 ± 7.8 ms	1 ± 0.85
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	2.93 ± 0.37 ms	2.88 ± 0.39 ms	1.02 ± 0.19
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	8.75 ± 1.5 ms	8.62 ± 1.1 ms	1.02 ± 0.22
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	1.54 ± 0.21 ms	1.53 ± 0.22 ms	1.01 ± 0.2
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.15 ± 3.3 ms	3.16 ± 3.4 ms	0.996 ± 1.5
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	5.51 ± 6.8 ms	5.31 ± 6.6 ms	1.04 ± 1.8
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.12 ± 0.033 s	0.117 ± 0.0067 s	1.03 ± 0.29
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.154 ± 0.11 s	0.147 ± 0.092 s	1.05 ± 0.97
Simple Pendulum/OOP/Shooting(Tsit5())	0.546 ± 0.31 ms	0.547 ± 0.35 ms	0.997 ± 0.86
time_to_load	4.47 ± 0.047 s	4.51 ± 0.026 s	0.992 ± 0.012

### 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).

[github-actions]

[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

	@views function __mirk_loss!(resid, u, p, y, pt::TwoPointBVProblem, bc!::Tuple{BC1, BC2},
	residual, mesh, cache, _, trait::DiffCacheNeeded) where {BC1, BC2}
	residual, mesh, cache, _, trait::DiffCacheNeeded, constraint::Val{false}) where {BC1, BC2}
	@views function __mirk_loss!(
	resid, u, p, y, pt::TwoPointBVProblem, bc!::Tuple{BC1, BC2}, residual, mesh,
	cache, _, trait::DiffCacheNeeded, constraint::Val{false}) where {BC1, BC2}