Open HildingElmqvist opened 5 years ago
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.
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
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.
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]
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?
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).
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.
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:
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?
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.
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:
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