Let me summarize the original problem of this thread because I am not quite sure if it has been resolved completely:
If the function f(u,p,t) in the ODEProblem uses besides some “tunable” (or “to-be-determined”) parameters (assembled into p) also some parameters that are not to be tuned/estimated such as the (mostly time-varying) control inputs, it is unclear how to estimate the former set of parameters following the procedure described in the manual for DiffEqParamEstim.jl. It was unclear to the author of the original post – and I think it has not yet been cleared out – how to tell the optimizer (to be called) over which parameters the optimization should be performed (and which parameters should be left intact).
It appears that the author @Anthony of the original post is happy with the advice to solve the problem by defining the control signal as a function of time. With this solution he soon runs into troubles with pulsed control signals, as he describes in his next post. Obviously the significant discontinuiuty (of the sharp edges) is the culprit. It seems to be possible to mittigate the problems by specifying the times of the switchnig through the tstops property (or perhaps even defining a callback that “does nothing”, if I understand it correctly). The problem seems to be answered/resolved.
Nonetheless, there may still be good reasons for defining the control input through parameters (possibly in combination with callback if the control evolves in time) – a feedback control is a prominent example. Or this is at least my understanding after following the discussion at DiffEqs: Hybrid (continuous+discrete) system / periodic callback - #38 by ChrisRackauckas. Is it possible to estimate just a subset of the parameters then?
The original post did not contain an example, that is why I am providing one. Perhaps it may be easier to continue a discussion then:
We consider a system modelled by \dot x(t) = a x(t) + b u(t) , where x is the state variable, u is the (known) control input. a and b are the parameters to be determined. (Note that I simplified the original problem by removing the offset parameter and considering just a perfectly linear system, for which alternative ways to solve the problem are available).
When generating the input-output data, I have also included some noise(s) in the model. First, a noise \zeta that affects the dynamics (aka random disturbance or process noise). And then some measurement noise \nu that adds to the measured state. Furthermore, the (noisy) state measurements are only available at some equidistant time instants with period T_s. Here is a full model used for generating the input-output data:
The actual control signal u will be designed in advance. It will consist of two (rectangular) pulses – first it goes from zero to one, then from one to minus one and then, finally, to zero again.
Below is a code for generating the data (advices on how to do this or that better are welcome but here it really just servers the purpose of creating the “experimental” data):
using DifferentialEquations
using Plots
# Setting up the continuous-time system
f(x,p,t) = p[1]*x + p[2]*p[3] # The standard RHS for an ODE, also the drift term for an SDE
g(x,u,p) = 0.1 # The diffusion term for an SDE
x₀ = 0.0 # The initial state
tspan = (0.0,10.0) # Time interval for the solution
a = -10.0 # Nominal parameters of the model
b = 10.0 # ---//---
u = 0.0 # The control is initially zero
p = [a,b,u]
prob = SDEProblem(f,g,x₀,tspan,p) # Formulating the problem as a stochastic ODE to include noisee
# Accounding for the external (control) input through callbacks
tₛ = [1.0,3.0,5.0] # Times when the control switches
condition1(u,t,integrator) = t==1.0 # The control switches to +1
affect1!(integrator) = integrator.p[3] = 1.0
cb1 = DiscreteCallback(condition1,affect1!)
condition2(u,t,integrator) = t==3.0 # The control switches to -1
affect2!(integrator) = integrator.p[3] = -1.0
cb2 = DiscreteCallback(condition2,affect2!)
condition3(u,t,integrator) = t==5.0 # The control switches back to zero again
affect3!(integrator) = integrator.p[3] = 0.0
cb3 = DiscreteCallback(condition3,affect3!)
cb = CallbackSet(cb1,cb2,cb3)
# Solving the problem
sol = solve(prob,SRIW1(),callback=cb,tstops=tₛ)
# Extracting the states at regular periods to pretend we measure them periodically
Tₛ = 0.2 # Sampling period
tₘ = range(tspan[1],tspan[2],step=Tₛ)
solₘ = sol(tₘ)
# Adding a "measurement" noise to sampled values
yₘ = solₘ.u + 0.1*randn(Float64,length(solₘ))
# (Preparing for) plotting the simulated states
p₁ = plot(sol, lw=2, label="Simulated in continuous time", xaxis = "t", yaxis = "x(t)")
scatter!(tₘ,yₘ,label="Noisy measurements in discrete time")
# Creating the sampled sequence of controls since these were not stored during simulation
function doublet(t)
if t < tₛ[1]
u = 0.0
elseif t < tₛ[2]
u = 1.0
elseif t < tₛ[3]
u = -1.0
else
u = 0.0
end
return u
end
uₘ = doublet.(tₘ) # Sampled controls
# (Preparing for) plotting the controls
tₚ = range(tspan[1],tspan[2],step=0.001) # Dense time for plotting the continuous-time controls
p₂ = plot(tₚ,doublet.(tₚ),xaxis="t",yaxis="u(t)", label="For simulation in continuous-time")
scatter!(tₘ,uₘ,label="Sampled in discrete-time")
# Plotting both the states and controls
plot(p₁,p₂,layout=(2,1))
Now, the input-output data set is ready: the vectors uₘ and yₘ (obtained with the sampling period Tₛ, possibly the vector tₘ of times when measurements were taken is available too). The parameters to be determined are a and b, and these will appear in the ODEProblem as p[1] and p[2] parameters.
prob4opt = ODEProblem(f,x₀,tspan,p)
But the optimization is certainly not to performed over p[3]. But how to include this information into the definition of the cost function?
cost_function = build_loss_objective(prob4opt,Tsit5(),L2Loss(tₘ,yₘ),
maxiters=10000,verbose=false)
Is it actually possible to make cost_function just a function of p[1] and p[2]? I was thinking about introducing equality constraints on p[3] in the optimization but that is nontrivial because the parameter (the control signal) evolves in time, it is not fixed.