# Issue with recreating MTK problems
**URL:** https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit
**Created:** [August 4, 2025, 5:23pm UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344 "2025-08-04T17:23:28Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@VPBML](https://discourse.julialang.org/u/VPBML)
#### Post date: [August 4, 2025, 5:23pm UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344/1 "2025-08-04T17:23:28Z")
Hi,
I am playing with a simple system consisting of a first-order dynamics and a PI-controller. The goal is to find a good tuning of the controller using optimization, with Integral Square Error as the objective function. I have made a code based on the [tutorial](https://docs.sciml.ai/ModelingToolkit/dev/examples/remake/).
I have observed that the BFGS solver always return the initial guess as the solution `sol1.u` with `sol1.retcode = 1` (even if the initial guess has wrong sign) and that the value of the objective function `sol1.objective` always has the same value although the solution has changed. The issue seems to be that although `newprob` gets the new parameters’ values from `x`, `sol = solve(newprob...)` still gives the same results as `sol0` (as if there is no change in parameters’ values).
I wonder what have I done wrong here and how to fix it. Thank you!
```julia-auto
using Optimization
using OptimizationOptimJL
using SymbolicIndexingInterface
using SymbolicIndexingInterface: parameter_values, state_values
using SciMLStructures: Tunable, canonicalize, replace, replace!
using PreallocationTools
using ModelingToolkit, Symbolics, OrdinaryDiffEq
import ModelingToolkit: t_nounits as t, D_nounits as D
@mtkmodel SimpleProcess begin
@components begin
process = FirstOrder(T = 5.0)
SP = Step(offset = 1.0, start_time = 5.0)
controller = PI(k = 1.0, T = 1.0, int.x = 0.0)
error = Add(k1 = 1.0, k2 = -1.0)
square_error = Power()
square_exponent = Constant(k = 2)
ISE = Integrator(k = 1.0, x = 0.0)
end
@equations begin
connect(SP.output, error.input1)
connect(process.output, error.input2)
connect(error.output, controller.err_input)
connect(controller.ctr_output, process.input)
connect(error.output, square_error.base)
connect(square_exponent.output, square_error.exponent)
connect(square_error.output, ISE.input)
end
end
odesys = SimpleProcess(; name = :odesys)
odesys = structural_simplify(odesys)
odeprob = ODEProblem(
odesys, [odesys.controller.int.x => 0.0, odesys.process.x => 0.0], (0.0, 50.0))
timesteps = 0.0:0.1:50.0
sol0 = solve(odeprob, Tsit5(); saveat = timesteps)
function loss(x, p)
odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
ps = parameter_values(odeprob) # obtain the parameter object from the problem
diffcache = p[4]
# get an appropriately typed preallocated buffer to store the `x` values in
buffer = get_tmp(diffcache, x)
# copy the current values to this buffer
copyto!(buffer, canonicalize(Tunable(), ps)[1])
# create a copy of the parameter object with the buffer
ps = replace(Tunable(), ps, buffer)
# set the updated values in the parameter object
setter = p[3]
setter(ps, x)
# remake the problem, passing in our new parameter object
newprob = remake(odeprob; p = ps)
timesteps = p[2]
sol = solve(newprob, AutoTsit5(Rosenbrock23()); saveat = timesteps)
# Handle solver failures
if !SciMLBase.successful_retcode(sol.retcode)
return 1e12 # Return a very high cost
end
data_vector = sol[odesys.ISE.output.u]
# Check for Inf or NaN values, which indicate an unstable solution.
if any(isinf, data_vector) || any(isnan, data_vector)
return 1e12 # Penalize unstable solutions heavily
end
return data_vector[end]
end
# manually create an OptimizationFunction to ensure usage of `ForwardDiff`, which will require changing the types of parameters from `Float64` to `ForwardDiff.Dual`
optfn = OptimizationFunction(loss, Optimization.AutoForwardDiff())
# function to set the parameters we are optimizing
setter = setp(odeprob, [odesys.controller.k, odesys.controller.T])
# `DiffCache` to avoid allocations. `copy` prevents the buffer stored by `DiffCache` from aliasing the one in `parameter_values(odeprob)`.
diffcache = DiffCache(copy(canonicalize(Tunable(), parameter_values(odeprob))[1]))
# parameter object is a tuple, to store differently typed objects together
optprob = OptimizationProblem(
optfn, [-1.0, 5.0], (odeprob, timesteps, setter, diffcache),
lb = [-Inf, 0.1], ub = [Inf, Inf])
sol1 = solve(optprob, BFGS())
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 4, 2025, 10:49pm UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344/2 "2025-08-04T22:49:53Z")
Which version is this?
---
### Author:  [@VPBML](https://discourse.julialang.org/u/VPBML)
#### Post date: [August 4, 2025, 11:09pm UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344/3 "2025-08-04T23:09:04Z")
Hi Chris, here are the versions of the packages:
- ModelingToolkit v9.76.0
- Symbolics v6.38.0
- OrdinaryDiffEq v6.95.1
- Optimization v4.5.0
- OptimizationOptimJL v0.4.3
- SymbolicIndexingInterface v0.3.42
- SciMLStructures v1.7.0
- PreallocationTools v0.4.29
My Julia version is 1.11.1. The example from [tutorial](https://docs.sciml.ai/ModelingToolkit/dev/examples/remake/) seems to work alright though.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 5, 2025, 1:17am UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344/4 "2025-08-05T01:17:11Z")
Open an issue, indeed it’s not entirely clear at first glance why this isn’t working like the tutorial
---
### Author:  [@VPBML](https://discourse.julialang.org/u/VPBML)
#### Post date: [August 5, 2025, 8:45am UTC](https://discourse.julialang.org/t/issue-with-recreating-mtk-problems/131344/5 "2025-08-05T08:45:00Z")
Thanks, Chris! I have opened an [issue](https://github.com/SciML/ModelingToolkit.jl/issues/3870).