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.
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!
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())