I’d like to optimize parameters in the MOL example shown here with an RMSD cost function to minimize.

Following a WIP example here to build the `OptimizationFunction`

, I have the following, but the results come back predicting the max of the bounded range.

```
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
using Plots
using DifferentialEquations, Optimization, OptimizationPolyalgorithms, Zygote, OptimizationOptimJL
using SciMLSensitivity
using Statistics
@parameters t x
p = @parameters Dn, Dp
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
eqs = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x),
Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
bcs = [u(0,x) ~ sin(pi*x/2),
v(0,x) ~ sin(pi*x/2),
u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0,1.0)]
@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2.0])
discretization = MOLFiniteDifference([x=>0.1],t)
pde_prob = discretize(pdesys,discretization)
sol = solve(pde_prob,Tsit5(),saveat=1.0)
meas = sol.u[2] .+ 0.01*randn(18)
idcs = Int.(ModelingToolkit.varmap_to_vars([p[1]=>1, p[2]=>2], p))
function loss(θ)
new_pde_prob = remake(pde_prob,p=θ[idcs])
sol = solve(new_pde_prob, Tsit5())
return sqrt(mean( (sol.u[2] .- meas).^2 )), Array(sol)
end
ad_type = Optimization.AutoZygote()
opt_func = OptimizationFunction((x,p)->loss(x), ad_type)
guess = [1.0,1.0]
opt_prob = OptimizationProblem(opt_func, guess, lb=[0.1,0.1],ub=[4.0,4.0])
res = Optimization.solve(opt_prob, GradientDescent())
```

Here’s the cost landscape in the bounded region:

And here’s the predicted parameters:

```
u: 2-element Vector{Float64}:
3.999999997189557
3.9999999981732874
```