Discuss using MTK in loss functions in the FAQ

Fixes #1692

Author: ChrisRackauckas

docs/src/basics/FAQ.md

@@ -41,5 +41,58 @@ ERROR: TypeError: non-boolean (Num) used in boolean context
 
 then it's likely you are trying to trace through a function which cannot be
 directly represented in Julia symbols. The techniques to handle this problem,
-such as `@register_symbolic`, are described in detail 
+such as `@register_symbolic`, are described in detail
 [in the Symbolics.jl documentation](https://symbolics.juliasymbolics.org/dev/manual/faq/#Transforming-my-function-to-a-symbolic-equation-has-failed.-What-do-I-do?-1).
+
+## Using ModelingToolkit with Optimization / Automatic Differentiation
+
+If you are using ModelingToolkit inside of a loss function and are having issues with
+mixing MTK with automatic differentiation, getting performance, etc... don't! Instead, use
+MTK outside of the loss function to generate the code, and then use the generated code
+inside of the loss function.
+
+For example, let's say you were building ODEProblems in the loss function like:
+
+```julia
+function loss(p)
+    prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```
+
+Since `ODEProblem` on a MTK `sys` will have to generate code, this will be slower than
+caching the generated code, and will required automatic differentiation to go through the
+code generation process itself. All of this is unnecessary. Instead, generate the problem
+once outside of the loss function, and remake the prob inside of the loss function:
+
+```julia
+prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+function loss(p)
+    remake(prob,p = ...)
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```
+
+Now, one has to be careful with `remake` to ensure that the parameters are in the right
+order. One can use the previously mentioned indexing functionality to generate index
+maps for reordering `p` like:
+
+```julia
+p = @parameters x y z
+idxs = ModelingToolkit.varmap_to_vars([p[1] => 1, p[2] => 2, p[3] => 3], p)
+p[idxs]
+```
+
+Using this, the fixed index map can be used in the loss function. This would look like:
+
+```julia
+prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+idxs = Int.(ModelingToolkit.varmap_to_vars([p1 => 1, p2 => 2], p))
+function loss(p)
+    remake(prob,p = p[idxs])
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```