Nonlinear Least Squares Problems with ModelingToolkit

We need to better document it. This works:

using NonlinearSolve
using ModelingToolkit

@variables begin
    cₜ
    k
    ϕₘ, [bounds = (1, Inf)]
    η, [bounds = (0, Inf)]
    b, [bounds = (0, Inf)]
end

@parameters begin
    t[1:100]
    ϕ[1:100]
end

krieger = @. (1 - ϕ/ϕₘ)^(-η * ϕₘ)
eqs = @. 0 ~ b + cₜ * t + k * krieger

@named ns = NonlinearSystem(eqs, [cₜ, k, ϕₘ, η, b], [t, ϕ])
ns = ModelingToolkit.complete(structural_simplify(ns, fully_determined = false))

prob = NonlinearLeastSquaresProblem(ns, [k => 0.0
                                ϕₘ => 0.0
                                η => 0.0
                                b => 0.0],
                                [
                                    t => rand(100)
                                    ϕ => rand(100)
                                ])
sol = solve(prob)

We probably want to flip that default on NonlinearSystem to allow for unbalanced. That is only supposed to be a check on ODESystem.