DifferentialEquations & inputs

General Usage
1

I’m trying to figure out the best way to pose a (control) system model of form:

\frac{dx}{dt} = f(x,u,t;\theta)

and use this model with DifferentialEquations. The “problem” is to merge two contradicting constraints:

  • DifferentialEquations requires the model to be posed as \frac{dx}{dt} = f(x,t;\theta), and
  • I would like to be able to pose the model in the form with input u, so that I can use Automatic Differentiation to linearize the model.

To test my idea, I used the “academic” model of level h [h] in a water tank, with input \dot{m}_\mathrm{i} [mdi], and effluent \dot{m}_\mathrm{e} driven by gravity, \dot{m}_\mathrm{e} = K\sqrt{h/h^\varsigma} where K is some valve constant, and h^\varsigma [hs] is some scaling level to keep the square root argument dimensionless. With density \rho [rho] and cross sectional area A [A], the model for the system is:

\frac{dh}{dt} = \frac{\dot{m}_\mathrm{i}-K\sqrt{h/h^\varsigma}}{\rho A}

I use the following parameters:

  • K=5, h^\varsigma = 3, \rho = 1, A=5.

I use the following initial value:

  • h(0) = 1.5.

At first, I want to set \dot{m}_\mathrm{i} = 2. Later on, I want to change this input value, e.g., by setting \dot{m}_\mathrm{i} = 2 + 0.1\sin(2\pi t/10), or by using a proportional controller \dot{m}_\mathrm{i} = K_\mathrm{p}(h_\mathrm{ref} - h). The point is: I want to make it as flexible as possible, while also allowing for Automatic Differentiation.

So… first I created a parameter tuple:

p = (rho=1.,A=5.,K=5.,hs=3.)

Next, I created a model function with input:

function mod_i(h,mdi,p,t)
    dh = (mdi-p.K*sqrt(h/p.hs))/(p.rho*p.A)
end

With an operating point (h^*=1.5,\dot{m}_\mathrm{i}^*=2), I can now linearize the model using ForwardDiff:

julia> ForwardDiff.derivative(h -> mod_i(h,2,p,0),1.5)
-0.2357022603955158
julia> ForwardDiff.derivative(mdi -> mod_i(1.5,mdi,p,0),2)
0.2

So far, so good. Next, I need to create a model function suitable for DifferentialEquations… and I attempt:

function mod(h,p,t)
    mdi = 2.
    return mod_i(h,mdi,p,t)
end

and set up the simulation case as:

h0 = 1.5
tspan = (0.,15.)
prob = ODEProblem(mod,h0,tspan)
sol = solve(prob)

Problem is: Julia crashes “big time”, with error message:

type Nothing has no field K
......

QUESTION:

  • What have I done wrong? Why does Julia crash?
1 Like
2

you just forgot to pass the parameters to the ODEProblem, like so

prob = ODEProblem(mod,h0,tspan,p)
2 Likes
3

Ooops… I did it again…

1 Like