Automatic differentiation of functions called within equations

[HildingElmqvist]

DAE index reduction requires that certain equations are differentiated versus time. The symbolic module of Modia contains partial differentiation of certain standard functions. A more general approach should be introduced for Modia. In particular, there will be many different fluid media functions calculating medium properties which needs to be differentiated.

The approach for Modelica is introduced in the paper Using Automatic Differentiation for Partial Derivatives of Functions in Modelica.

A good little first example might be the straightFlanksCam function which I translated to Julia:

function straightFlanksCam(
    theta::Float64;
    R1::Float64=1.0,   # Base circle radii
    R2::Float64=0.5,   # Nose radii
    d::Float64=2.0     # Centre distance
    )
    thetamod = atan(sin(theta), cos(theta)) # to get angle in interval -pi..pi
    fi0 = asin((R1 - R2)/d)
    if thetamod > 0
        fi1 = pi/2 - fi0 - thetamod
    else
        fi1 = - pi/2 + fi0 - thetamod
    end 
    fi2 = atan(R2*cos(fi0), d + R2*sin(fi0))
    if abs(thetamod) > pi/2-fi0
        x = R1*cos(theta)
        y = R1*sin(theta)
    elseif abs(thetamod) >= fi2
        L = R1/cos(fi1)
        x = L*cos(theta)
        y = L*sin(theta)
    else
        L = d*cos(abs(thetamod)) + sqrt(R2^2 - d^2*sin(abs(thetamod))^2)
        x = L*cos(theta)
        y = L*sin(theta)
    end
    return r = [x, y]
end

This function is not used for index reduction but shows another need, i.e. finding the direction of the force on a body, by finding the normal to a parametric surface.

Using AD techniques could initially be tested without integration in Modia.

@MartinOtter will provide other examples involving media functions.

CC: @crlaugh, @jrevels

[HildingElmqvist]

It is also of interest to study the Modelica specification

The section: 12.7 Declaring Derivatives of Functions (page 160) describes how one can provide derivative functions manually. It shows several examples of what the interfaces to derivative functions would be.

[HildingElmqvist]

The current Modia implementation regarding differentation of non-standard functions uses a naming convention.

The tests for symbolic transformation illustrates the current convention:

  der = showDifferentiate(:(y = f(x, 5, z)))
  @test der == "der(y) = f_der_1(x, 5, z) * der(x) + f_der_3(x, 5, z) * der(z)"
  
  der = showDifferentiate(:(y = f(x, 5, g(z))))
  @test der == "der(y) = f_der_1(x, 5, g(z)) * der(x) + f_der_3(x, 5, g(z)) * (g_der(z) * der(z))"

The code for this handling is here:

            # der(f(e1, e2, e3)) = f_der_1(e1, e2, e3)*der(e1) + f_der_2(e1, e2, e3)*der(e2) + f_der_3(e1, e2, e3)*der(e3)
            arguments = e.args[2:end]
            diff = Expr(:call, +)
            
            for i in 1:length(arguments)
                if length(arguments) == 1
                    f_der = Symbol(string(op, "_der"))
                else 
                    f_der = Symbol(string(op, "_der_", i))
                end
            
                f_der_i = Expr(:call, f_der)
                for a in arguments
                    push!(f_der_i.args, a)
                end
            
                term = mult(f_der_i, differentiate(arguments[i])) 
                if term != zero
                    push!(diff.args, term)      
                end
            end       

[jrevels]

An example calculating the Jacobian of the OP using ForwardDiff + StaticArrays:

# only two changes:
# 1. remove argument type constraints
# 2. return a tuple instead of an array (this isn't absolutely necessary, but is more efficient)
function straightFlanksCam(theta; R1 = 1.0, #= Base circle radii =#
                           R2 = 0.5, #= Nose radii =#
                           d = 2.0 #= Centre distance =#)
    thetamod = atan(sin(theta), cos(theta)) # to get angle in interval -pi..pi
    fi0 = asin((R1 - R2)/d)
    if thetamod > 0
        fi1 = pi/2 - fi0 - thetamod
    else
        fi1 = - pi/2 + fi0 - thetamod
    end 
    fi2 = atan(R2*cos(fi0), d + R2*sin(fi0))
    if abs(thetamod) > pi/2-fi0
        x = R1*cos(theta)
        y = R1*sin(theta)
    elseif abs(thetamod) >= fi2
        L = R1/cos(fi1)
        x = L*cos(theta)
        y = L*sin(theta)
    else
        L = d*cos(abs(thetamod)) + sqrt(R2^2 - d^2*sin(abs(thetamod))^2)
        x = L*cos(theta)
        y = L*sin(theta)
    end
    return x, y
end

using ForwardDiff, StaticArrays

function straightFlanksCamJacobian(theta; R1 = 1.0, R2 = 0.5, d = 2.0)
    x = SVector(theta, R1, R2, d)
    f = x -> SVector(straightFlanksCam(x[1], R1=x[2], R2=x[3], d=x[4]))
    return ForwardDiff.jacobian(f, x)
end

Then, for example:

julia> straightFlanksCamJacobian(rand())
2×4 SArray{Tuple{2,4},Float64,2,8}:
 -0.25304  0.0  1.00308    0.997953
  2.49335  0.0  0.0202763  0.0201727

If you have an analytical derivative it'd be good to test against that, or we could just test against finite differencing.

[HildingElmqvist]

Excellent: Some test cases:

straightFlanksCamJacobian(0.0)*[1; 0; 0; 0] == [0.0, 2.5]
straightFlanksCamJacobian(pi/2)*[1; 0; 0; 0] == [-1.0, 6.123233995736766e-17]
straightFlanksCamJacobian(pi)*[1; 0; 0; 0] == [-1.2246467991473532e-16, -1.0]
straightFlanksCamJacobian(3*pi/2)*[1; 0; 0; 0] == [ 1.0, -1.8369701987210297e-16]

[HildingElmqvist]

It seems the user function needs to be wrapped in order to have a function with a single vector input for ForwardDiff. What would be the best way of automatically provide the wrapping function (without eval())? Or could a generic wrapping function be constructed?

[MartinOtter]

Here is a very simple model that requires differentiation of medium function specificEnthalpy_water:

module Test_MediaDifferentiation1
    using  Modia

    const cp_const = 4184
    const T0       = 273.15

    specificEnthalpy_water(T) = cp_const*(T - T0)

    @model Test_Differentiation1 begin 
        T  = Float(start=300.0, state=true)
        h  = Float(start=0.0, state=false)

        @equations begin
            h = specificEnthalpy_water(T)
            der(h) + der(T) = time 
        end
    end

    result = simulate(Test_Differentiation1, 0.1)
end

Simulating model: Test_Differentiation1
ERROR: LoadError: MethodError: Cannot `convert` an object of type Nothing to an object of type Int64
Closest candidates are:
  convert(::Type{T<:Real}, ::Unitful.Quantity) where T<:Real at D:\otter\home\.julia\packages\Unitful\MlDo0\src\conversion.jl:138
  convert(::Type{T<:Real}, ::Unitful.Level) where T<:Real at D:\otter\home\.julia\packages\Unitful\MlDo0\src\logarithm.jl:36
  convert(::Type{T<:Real}, ::Unitful.Gain) where T<:Real at D:\otter\home\.julia\packages\Unitful\MlDo0\src\logarithm.jl:94
  ...
Stacktrace:
 [1] setindex!(::Array{Int64,1}, ::Nothing, ::Int64) at .\array.jl:769
 [2] copyto! at .\abstractarray.jl:731 [inlined]
 [3] copyto! at .\abstractarray.jl:723 [inlined]
 [4] Type at .\array.jl:497 [inlined]
 [5] convert at .\array.jl:489 [inlined]
 [6] push!(::Array{Array{Int64,1},1}, ::Array{Nothing,1}) at .\array.jl:855
 [7] handleSingularities(::Array{Any,1}, ::Array{Any,1}, ::Dict{Symbol,Int64}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Int64,1}, ::Array
{Any,1}, ::Array{Any,1}, ::Bool) at Modia\src\symbolic\DAEquations\BasicStructuralTransform.jl:521

Running this with Modia, gives an error (as to be expected; it is just a bit strange that the error occurs in Unitful).

I guess Hilding (@HildingElmqvist) you need to have a look, before passing the differentiation problem to Jarrett (@jrevels).

[MartinOtter]

I have updated the example, so that it includes the derivative of the medium function with ForwardDiff using the naming convention sketched by Hilding above:

module Test_MediaDifferentiation1
    using  Modia
    using ForwardDiff 

    const cp_const = 4184
    const T0       = 273.15

    specificEnthalpy_water(T) = cp_const*(T - T0)

    specificEnthalpy_water_der_1(T,T_der_1) = 
            ForwardDiff.derivative(T->specificEnthalpy_water(T),T)*T_der_1

    @model Test_Differentiation1 begin 
        T  = Float(start=300.0, state=true)
        h  = Float(start=0.0, state=false)

        @equations begin
            h = specificEnthalpy_water(T)
            der(h) + der(T) = time 
        end
    end

    result = simulate(Test_Differentiation1, 0.1)
end

Simulating model: Test_Differentiation1
ERROR: LoadError: MethodError: Cannot `convert` an object of type Nothing to an object of type Int64

This gives still the error from above. Hilding (@HildingElmqvist) can you have a look.

[MartinOtter]

I have experimented a bit with a tiny set of media functions, computing their derivatives and plotting them (see below). It seems we should implement a "ModiaMedia" package that is completely independent of Modia and just provides the media models in form of Julia functions (along the sketch of Test_Media.jl). The goal would be to convert Modelica.Media to Julia in a semi-automatic way and make improvements on the way (e.g. store the Modelica.Media.IdealGases.Common.SingleGasesData not as constants but as json file, that is only read, when gas data from this file is needed).

module Test_MediaDifferentiation2

using  StaticArrays
import ModiaMath
using  ModiaMath.Unitful

### Data structure for ideal gases

struct DataRecord
    name::AbstractString       # Name of ideal gas
    MM::Float64                # SI.MolarMass "Molar mass"
    Hf::Float64                # SI.SpecificEnthalpy "Enthalpy of formation at 298.15K"
    H0::Float64                # SI.SpecificEnthalpy "H0(298.15K) - H0(0K)"
    Tlimit::Float64            # SI.Temperature Tlimit "Temperature limit between low and high data sets"
    alow::SVector{7,Float64}   # "Low temperature coefficients a"
    blow::SVector{2,Float64}   # "Low temperature constants b"
    ahigh::SVector{7,Float64}  # "High temperature coefficients a"
    bhigh::SVector{2,Float64}  # "High temperature constants b"
    R::Float64                 # SI.SpecificHeatCapacity R "Gas constant"

    DataRecord(;name=nothing,
                MM=nothing,
                Hf=nothing,
                H0=nothing,
                Tlimit=nothing,
                alow=nothing,
                blow=nothing,
                ahigh=nothing,
                bhigh=hothing,
                R=nothing) = 
                new(name,MM,Hf,H0,Tlimit,alow,blow,ahigh,bhigh,R)
end



### Enumerations and constants

"""
    @enum referenceChoice

Enumeration defining the reference enthalpy of a medium. Possible values:

- `ZeroAt0K`: The enthalpy is 0 at 0 K (default), if the enthalpy of formation is excluded
- `ZeroAt25C`: The enthalpy is 0 at 25 degC, if the enthalpy of formation is excluded
- `UserDefined`: The user-defined reference enthalpy is used at 293.15 K (25 degC)
"""
@enum ReferenceChoice ZeroAt0K  ZeroAt25C  UserDefined

const excludeEnthalpyOfFormation = true       # If true, enthalpy of formation Hf is not included in specific enthalpy hc
const referenceChoice            = ZeroAt0K   # Choice of reference enthalpy
const h_offset                   = 0.0        # User defined offset for reference enthalpy, if referenceChoice = UserDefined



### Data records for specific ideal gases

const H2O = DataRecord(
    name="H2O",
    MM=0.01801528,
    Hf=-13423382.81725291,
    H0=549760.6476280135,
    Tlimit=1000,
    alow=SVector(-39479.6083,575.573102,0.931782653,0.00722271286,-7.34255737e-006,
        4.95504349e-009,-1.336933246e-012),
    blow=SVector(-33039.7431,17.24205775),
    ahigh=SVector(1034972.096,-2412.698562,4.64611078,0.002291998307,-6.836830479999999e-007,
        9.426468930000001e-011,-4.82238053e-015),
    bhigh=SVector(-13842.86509,-7.97814851),
    R=461.5233290850878)


const N2 = DataRecord(
    name="N2",
    MM=0.0280134,
    Hf=0,
    H0=309498.4543111511,
    Tlimit=1000,
    alow=SVector(22103.71497,-381.846182,6.08273836,-0.00853091441,1.384646189e-005,-9.62579362e-009,
        2.519705809e-012),
    blow=SVector(710.846086,-10.76003744),
    ahigh=SVector(587712.406,-2239.249073,6.06694922,-0.00061396855,1.491806679e-007,-1.923105485e-011,
        1.061954386e-015),
    bhigh=SVector(12832.10415,-15.86640027),
    R=296.8033869505308);




### Functions for ideal gases
h_T(data::DataRecord, T; 
    exclEnthForm::Bool         = excludeEnthalpyOfFormation,  # If true, enthalpy of formation Hf is not included in specific enthalpy h
    refChoice::ReferenceChoice = referenceChoice,             # Choice of reference enthalpy
    h_off::Float64             = h_offset                     # User defined offset for reference enthalpy, if referenceChoice = UserDefined
   ) = (T < data.Tlimit ?
          data.R*((-data.alow[1] + T*(data.blow[1] + data.alow[2]*log(T) + T*(
                    data.alow[3] + T*(0.5*data.alow[4] + T*(1/3*data.alow[5] + T*(0.25*data.alow[6] + 0.2*data.alow[7]*T))))))/T) :
          data.R*((-data.ahigh[1] + T*(data.bhigh[1] + data.ahigh[2]*log(T) + T*(
                    data.ahigh[3] + T*(0.5*data.ahigh[4] + T*(1/3*data.ahigh[5] + T*(0.25*data.ahigh[6] + 0.2*data.ahigh[7]*T))))))/T)) + 
         (exclEnthForm ? -data.Hf : 0.0) + 
         (refChoice == ZeroAt0K ? data.H0 : 0.0) + 
         (refChoice == UserDefined ? h_off : 0.0)

specificEnthalpy(data::DataRecord, T) = h_T(data,T)


### Determine derivatives

using ForwardDiff

specificHeatCapacityCp(data::DataRecord, T) = 
   ForwardDiff.derivative(T->specificEnthalpy(data,T),T)

specificEnthalpy_der_1(data::DataRecord, T, T_der_1) = 
   specificHeatCapacityCp(data::DataRecord, T)*T_der_1


### Plot medium data

T      = collect(300.0:10.0:2000.0)
h      = zeros(length(T))
cp     = zeros(length(T))
der_h  = zeros(length(T))


# Select gas
data = N2

for i in 1:length(T)
    h[i]      = specificEnthalpy(data,T[i])
    cp[i]     = specificHeatCapacityCp(data,T[i])
    
    T_der_1   = 10.0*i/length(T)
    der_h[i]  = specificEnthalpy_der_1(data,T[i],T_der_1)
end
mediumData = Dict{AbstractString,Any}()
mediumData["T"]         = T*1u"K"
mediumData["h"]         = h*1u"J/kg"
mediumData["cp"]        = cp*1u"J/(kg*K)"
mediumData["der(h)"]    = der_h*1u"J/(kg*s)"

ModiaMath.plot(mediumData, "h"     , xAxis="T", heading=data.name, figure=1)
ModiaMath.plot(mediumData, "cp"    , xAxis="T", heading=data.name, figure=2)
ModiaMath.plot(mediumData, "der(h)", xAxis="T", heading=data.name, figure=3)
end

One of the figures generated by the code above:

[HildingElmqvist]

There is an issue how to make the lookup of the derivative function.

Presently a new Symbol is created for the derivative function:
f_der = Symbol(string(op, "_der"))
and a new Expr is created:
f_der_i = Expr(:call, f_der)
and the same arguments are provided for the derivative function:

 for a in arguments
    push!(f_der_i.args, a)
end

However, the Symbol needs to be looked up in the user module if the user has to make the wrapper.

How can this look up be made in the correct module?

[HildingElmqvist]

The solution was to use the getfield() function. So now, using Modia master, derivative functions are found and called. But a wrapper needs to be made maually using, for example, the ForwardDiff package.