ModellingToolkit ODESystem definition throwing error

Specific Domains Numerics
,
1

I am defining and ODESystem using ModellingToolkit and I am running into a problem in the definition of the equations.

The code below is the function to define the ODESystem

function connectFFFB_Grad(plant, ref, pid, mlp,lookahead_times,N,plantMatrices;name=:FF_FB_G)   
    @parameters t 
    dt = Differential(t)
    @variables intError(t)
    # @variables dydw[1:N](t) #grad output wrt weights
    @variables dx1dw[1:N](t) #grad 4th plant state wrt weights
    @variables dx2dw[1:N](t) #grad 4th plant state wrt weights
    @variables dx3dw[1:N](t) #grad 4th plant state wrt weights
    @variables dx4dw[1:N](t) #grad 4th plant state wrt weights
    A,B,C,D = plantMatrices
    # @variables e(t)
    eqs = vcat(
            get_input(pid,1) ~ get_output(ref, 1)  - get_output(plant,1),
            get_output(pid,1) ~ get_input(plant, 1),
            get_output.((ref,), 2:length(lookahead_times)+1) .~ get_inputs(mlp,length(lookahead_times)),
            get_output(mlp,1) ~ get_input(plant,2),
            dt(intError) ~  abs2(get_output(ref,1)  - get_output(plant,1)),

            dt.(dx1dw) .~ A[1,:]*vcat(dx1dw,dx2dw,dx3dw,dx4dw), 
            dt.(dx2dw) .~ A[2,:]*vcat(dx1dw,dx2dw,dx3dw,dx4dw),
            dt.(dx3dw) .~ A[3,:]*vcat(dx1dw,dx2dw,dx3dw,dx4dw), 
            dt.(dx4dw) .~ A[4,:]*vcat(dx1dw,dx2dw,dx3dw,dx4dw), 
            # dt.(dydw) .~ C*[A*vcat(dx1dw,dx2dw,dx3dw,dx4dw).+B[:,2]*get_hiddens(mlp,N)], 
            # # e ~ abs2(get_output(ref,1)  - get_output(plant,1)),
    ) 
    return ODESystem(eqs; systems=[plant,ref,pid,mlp],name=name)
end

I am getting the following error

ERROR: MethodError: no method matching *(::Vector{Int64}, ::Vector{Num})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
  *(::StridedMatrix{T}, ::StridedVector{S}) where {T<:Union{Float32, Float64, ComplexF32, ComplexF64}, S<:Real} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:44
  *(::StridedVecOrMat{T} where T, ::LinearAlgebra.LQPackedQ) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lq.jl:278
  ...
Stacktrace:
 [1] connectFFFB_Grad(plant::ODESystem, ref::ODESystem, pid::ODESystem, mlp::ODESystem, lookahead_times::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, N::Int64, plantMatrices::NTuple{4, Matrix{Int64}}; name::Symbol)

I know that the error is in the te last four equations as I had previously function defining the same system excluding those four variables and it worked well.

It seems that the multiplication between the matrix column and the vector with the variables is not defined but I don’t understand why. Especially as I have another ODESystem defined by the following function that doesn’t pose any problems.

function linear_plant(A,B,C,D; name=:lp)
    
    num_inputs = size(B,2)
    num_states = size(A,1)
    num_outputs = size(C,1)

    @parameters t
    dt = Differential(t)

    @variables x[1:num_states](t)
    @variables u[1:num_inputs](t)
    @variables o[1:num_outputs](t)
    
    eqs = vcat(
                dt.(x) .~ A*x .+ B*u,
                o .~  C*x .+ D*u
    )    
    return ODESystem(eqs;name=name)   
end

A is a 4x4 matrix defined as follows

A = [0 1 0 0; -2 -2 1 0; 0 0 0 1; 0 0 -2 -2];

I would appreciate any help understanding this error.

2

I think you need a broadcast in there:

A[1,:] .* vcat(dx1dw,dx2dw,dx3dw,dx4dw)
3

Thanks for the suggestion contradict. I tried that and I get an error

julia> ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 4 and 40")
Stacktrace:
  [1] _bcs1
    @ ./broadcast.jl:501 [inlined]
  [2] _bcs
    @ ./broadcast.jl:495 [inlined]
  [3] broadcast_shape
    @ ./broadcast.jl:489 [inlined]
  [4] combine_axes
    @ ./broadcast.jl:484 [inlined]
  [5] _axes
    @ ./broadcast.jl:209 [inlined]
  [6] axes
    @ ./broadcast.jl:207 [inlined]
  [7] combine_axes
    @ ./broadcast.jl:484 [inlined]
  [8] instantiate
    @ ./broadcast.jl:266 [inlined]
  [9] materialize
    @ ./broadcast.jl:883 [inlined]
 [10] connectFFFB_Grad(plant::ODESystem, ref::ODESystem, pid::ODESystem, mlp::ODESystem, lookahead_times::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, N::Int64, plantMatrices::NTuple{4, Matrix{Int64}}; name::Symbol)

I think A[1,:] is a vector of length 4 and vcat(dx1dw,dx2dw,dx3dw,dx4dw) is 4xN in this case N=10.

4

Oh, I see, I misunderstood. It looks like there is no vector-matrix product defined for the Num type, I think the proposed solution for this is support for Symbolic Arrays: https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/123

I suppose you could write it out manually, but I don’t know of a nice solution to this.

5

This is all undergoing major changes with symbolic arrays:

https://symbolics.juliasymbolics.org/dev/manual/arrays/

merged like yesterday, and it’s getting to ModelingToolkit:

https://github.com/SciML/ModelingToolkit.jl/pull/989

This kind of thing will make it so you can say “this is an array operation” as a single symbolic command, and should help in these cases. I’ll just say, let’s get the major thrust of this work done and then open an issue if you still have some problems are moving to the new form.

6

Thanks Chris!