search

ModiaSim / Modia.jl

Modeling and simulation of multidomain engineering systems
MIT License
319 stars 37 forks source link

linearize!(model, analytic=true, ... #122

Open johhell opened 2 years ago

johhell commented 2 years ago

I tried to linearize a simple electrical system and received an error message when using option analytic=true

The problem occurs when getDerivatives! is called from ForwardDiff.Jacobian, then _leq_mode.residual_value is set to:

Any[Dual{ForwardDiff.Tag{ModiaLang.var"#modelToLinearize!#54"{SimulationModel{Float64,OrderedCollections.OrderedDict{Symbol,Any},OrderedCollections.OrderedDict{Symbol,Any},Float64}},Float64}}(-400001.0,-1000.0,-0.0,-400000.0,-0.0)]

which can not be handled by function LinearEquationsIteration.

generated function getDerivatives!

code = quote
    function getDerivatives(_der_x, _x, _m, _time)::Nothing
        _m.time = ModiaLang.getValue(_time)
        _m.nGetDerivatives += 1
        instantiatedModel = _m
        _p = _m.evaluatedParameters
        _leq_mode = nothing
        time = _time
        var"T1.X.ixRe" = _x[1]
        var"T1.X.ixIm" = _x[2]
        var"C1.voltsRe" = _x[3]
        var"C1.voltsIm" = _x[4]
        var"T1.Omegarated" = _p[:Omegarated]
        var"T1.wReference" = _p[:wReference]
        begin
            local var"T1.X.voltsRe", var"T1.R.voltsRe", var"C1.ampsRe", var"T1.X.iparallelRe"
            _leq_mode = _m.linearEquations[1]
            _leq_mode.mode = -2
             ModiaBase.TimerOutputs.@timeit _m.timer "LinearEquationsIteration" 
                while ModiaBase.LinearEquationsIteration(_leq_mode, _m.isInitial, _m.time, _m.timer)
                    var"T1.X.voltsRe" = _leq_mode.vTear_value[1]
                    var"T1.R.voltsRe" = -1var"C1.voltsRe" + -1var"T1.X.voltsRe"
                    var"C1.ampsRe" = var"T1.R.voltsRe" / ((_p[:T1])[:R])[:r]
                    var"T1.X.iparallelRe" = -((var"T1.X.ixRe" + -1var"C1.ampsRe"))
                    _leq_mode.residual_value[1] = ustrip(((_p[:T1])[:X])[:rparallel] * var"T1.X.iparallelRe" - var"T1.X.voltsRe")
                end
            _leq_mode = nothing
        end
        var"der(T1.X.ixRe)" = -((var"T1.X.voltsRe" - ((_p[:T1])[:X])[:x] * -(var"T1.wReference" * var"T1.X.ixIm"))) / -(((_p[:T1])[:X])[:x] * (1 / var"T1.Omegarated"))
        begin
            local var"T1.X.voltsIm", var"T1.R.voltsIm", var"C1.ampsIm", var"T1.X.iparallelIm"
            _leq_mode = _m.linearEquations[2]
            _leq_mode.mode = -2
             ModiaBase.TimerOutputs.@timeit _m.timer "LinearEquationsIteration" while ModiaBase.LinearEquationsIteration(_leq_mode, _m.isInitial, _m.time, _m.timer)
                    var"T1.X.voltsIm" = _leq_mode.vTear_value[1]
                    var"T1.R.voltsIm" = -1var"C1.voltsIm" + -1var"T1.X.voltsIm"
                    var"C1.ampsIm" = var"T1.R.voltsIm" / ((_p[:T1])[:R])[:r]
                    var"T1.X.iparallelIm" = -((var"T1.X.ixIm" + -1var"C1.ampsIm"))
                    _leq_mode.residual_value[1] = ustrip(((_p[:T1])[:X])[:rparallel] * var"T1.X.iparallelIm" - var"T1.X.voltsIm")
                end
            _leq_mode = nothing
        end
        var"der(T1.X.ixIm)" = -((var"T1.X.voltsIm" - ((_p[:T1])[:X])[:x] * (var"T1.wReference" * var"T1.X.ixRe"))) / -(((_p[:T1])[:X])[:x] * (1 / var"T1.Omegarated"))
        var"T1.R.p.vRe" = -((-1var"T1.R.voltsRe" + -1var"T1.X.voltsRe"))
        var"T1.R.p.vIm" = var"T1.R.voltsIm" + var"T1.X.voltsIm"
        var"C1.Omegarated" = _p[:Omegarated]
        var"C1.wReference" = _p[:wReference]
        var"der(C1.voltsRe)" = -((var"C1.ampsRe" * (_p[:C1])[:b] - -(var"C1.wReference" * var"C1.voltsIm"))) / -(1 / var"C1.Omegarated")
        var"der(C1.voltsIm)" = -((var"C1.ampsIm" * (_p[:C1])[:b] - var"C1.wReference" * var"C1.voltsRe")) / -(1 / var"C1.Omegarated")
        _der_x[1] = var"der(T1.X.ixRe)"
        _der_x[2] = var"der(T1.X.ixIm)"
        _der_x[3] = var"der(C1.voltsRe)"
        _der_x[4] = var"der(C1.voltsIm)"
        if _m.storeResult
            ModiaLang.addToResult!(_m, _der_x, time, var"T1.Omegarated", var"T1.wReference", var"T1.X.voltsRe", var"T1.X.voltsIm", var"T1.X.iparallelRe", var"T1.X.iparallelIm", var"T1.R.voltsRe", var"T1.R.p.vRe", var"T1.R.voltsIm", var"T1.R.p.vIm", var"C1.Omegarated", var"C1.wReference", var"C1.ampsRe", var"C1.ampsIm")
        end
        return nothing
    end
end

versions

"Modia" version = "0.5.0" "ModiaBase" version = "0.7.5" "ModiaLang" version = "0.8.1" Julia: Version 1.5.3 (2020-11-09) platform: LINUX x86_64

PS:

After reading your paper from the ongoing MODELICA conference, I had a better understanding but couldn't resolve it.

MartinOtter commented 2 years ago

Yes, I know. There are also other situations where analytic=true does not work. I spent some time on it, but did not find a solution. I have documented this in the description of linearize (analytic=true might not work) and used analytic=false (numeric linearization) as a default - which should always work. A workaround is to perform linearization numerically with Double64, which should give an even higher precision as analytic linearization with Float64.

model2 = SimulationModel{Double64}(model1)   # transform instantiated model1 from Float64 to Double64
(A,x0) = linearize!(model2)
johhell commented 2 years ago

below my proposal for linearize!(model, analytic=true) . Not perfect, but it works

Modification of file ModiaBase...EquationAndStateInfo.jl

structure: mutable struct LinearEquations{FloatType <: Real}

For Vectors/Matrix: A,x,b,residuals,luA , I changed the type from FloatType to Union{FloatType,Any} 

mutable struct LinearEquations{FloatType <: Real}
    odeMode::Bool                   # Set from the calling function after LinearEquations was instantiated (default: true)
                                    # = true (standard mode): Compute "x" from equation "residuals = A*x - b"
                                    # = false (DAE solver): During events (including initialization)
                                    #   compute "x" as for odeMode=true. Outside of events:
                                    #   (1) "x" is set from the outside (= der(x_dae) provided by DAE solver)
                                    #   (2) Compute "residuals" from "residuals := A*x - b"
                                    #   (3) From the outside copy "residuals" into "residuals_dae" of the DAE solver.

    A_is_constant::Bool             # = true, if A-matrix is constant
    x_names::Vector{String}         # Names of the x-variables
    x_lengths::Vector{Int}          # Lengths of the x-variables (sum(x_lengths) = length(x))
    x::Vector{Union{FloatType,Any}}            # Values of iteration variables
    nResiduals::Int                 # Number of residual variables

    A::Matrix{Union{FloatType,Any}}
    b::Vector{Union{FloatType,Any}}
    pivots::Vector{Int}             # Pivot vector if recursiveFactorization = true
    residuals::Vector{Union{FloatType,Any}}    # Values of the residuals FloatType vector; length(residuals) = sum(residuals_length) = sum(x_lengths)
    residual_value::AbstractVector  # Values of the residual variables ::Vector{Any}, length(residual_values) = nResiduals

    # Iteration status of for-loop
    mode::Int       # Operational mode (see function LinearEquationsIteration)
    niter::Int      # Current number of iterations in the fix point iteration scheme
    niter_max::Int  # Maximum number of iterations
    success::Bool   # = true, if solution of A*x = b is successfully computed
                    # = false, if solution is not computed; continue with fix point iteration

    # For warning message if niter > niter_max
    inconsistentPositive::Vector{String}
    inconsistentNegative::Vector{String}

    # Constructed during initialization
    residual_unitRanges::Vector{UnitRange{Int}} # residuals[residual_unitRanges[i]] = residual_value[i], if residual is a vector
    residual_indices::Vector{Int}               # residuals[residual_indices[i]] = residual_value[i], if residual is a scalar
    recursiveFactorization::Bool                # = true, if RecursiveFactorization.jl shall be used to solve the linear equation system

    luA::LU{Union{FloatType,Any},Array{Union{FloatType,Any},2}}       # lu-Decomposition of A
...

This works for nx==1, but not for nx>1 when lu! is called.

additional modification in RecursiveFactorization...lu.jl

zero() function for type Dual does not exist. So I changed to convert()

function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where Pivot
    m, n = size(A)
    minmn = length(ipiv)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
#   =====>>>>          amax = abs(zero(eltype(A)))   <<<<====
              amax = convert(eltype(A),0)
              for i = k:m
                  absi = abs(A[i,k])
                  if absi > amax
                      kp = i
...

version

name = "ModiaBase" version = "0.8.0-dev"