SciML / ModelingToolkitStandardLibrary.jl

A standard library of components to model the world and beyond
https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/
MIT License
119 stars 37 forks source link

`MatrixGain` array parameter does not work #213

Closed baggepinnen closed 1 year ago

baggepinnen commented 1 year ago

The example below used to work fine (MTKstdlib v1), but now errors with

julia> prob = ODEProblem(closed_loop, Pair[], (0., 10.))
ERROR: MethodError: Cannot `convert` an object of type 
  Bool to an object of type 
  Union{Float64, Int64, Matrix{Float64}}
using ModelingToolkit, LinearAlgebra
using ModelingToolkitStandardLibrary.Mechanical.Rotational
import ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: MatrixGain, Add, Step, Sine
t = Blocks.t
m1 = 1
m2 = 1
k = 1000 # Spring stiffness
c = 10   # Damping coefficient

@named inertia1 = Inertia(; J = m1)
@named inertia2 = Inertia(; J = m2)
@named spring = Spring(; c = k)
@named damper = Damper(; d = c)
@named torque = Torque(use_support=false)

function SystemModel(u=nothing; name=:model)
    eqs = [
        connect(torque.flange, inertia1.flange_a)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
    ]
    if u !== nothing 
        push!(eqs, connect(torque.tau, u.output))
        return @named model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper, u])
    end
    ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
end

@named d = Step(start_time=1, duration=2) # Disturbance
model = SystemModel() # Model with load disturbance
outputs = [inertia1.w, inertia2.w, inertia1.phi, inertia2.phi]
inputs = [torque.tau.u]
L = randn(1, 4)
@named state_feedback = MatrixGain(L)
@named add = Add() # To add the control signal and the disturbance

model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi]
connections = [
    [state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
    connect(add.input1, d.output)
    connect(add.input2, state_feedback.output)
    connect(add.output, model.torque.tau)
]
closed_loop = ODESystem(connections, t, systems=[model, state_feedback, add, d], name=:closed_loop)
closed_loop = structural_simplify(closed_loop)
prob = ModelingToolkit.ODEProblem(closed_loop, Pair[], (0., 10.))

The problem is related to the parameter MatrixGain.K, which in MTKstdlib v1 was a hard-coded constant, while in v2 it's a parameter

baggepinnen commented 1 year ago

@ven-k is it possible to define components with @mtkmodel and make the constructor take input arguments without making those arguments a parameter of the system?

baggepinnen commented 1 year ago

I can confirm that the error goes away with the old implementation of MatrixGain

ven-k commented 1 year ago

RN, any input to the component will have to be either a parameter or a variable. We can have constructor take in an input, here K, and get the nout = size(K, 2) and nin = size(K, 1) and pass them to the component. But nout and nin will be the parameters of the component.

baggepinnen commented 1 year ago

But nout and nin will be the parameters of the component.

This would be confusing since they cannot be modified the same way as standard parameters. In general, it would be good to let the @mtkmodel macro accept arguments like standard functions, the multibody library uses this a lot to select things like how to initialize variables, which orientation representation to use etc.

ven-k commented 1 year ago

The issue isn't because of MatrixGain; but it triggers it.

The error stems from the fact that the vs has array and hence the C is promoted to Union{Float64, Int64, Matrix{Float64}}. RN it doesn't handle when Boolean values are present.

In the above example, use_support from Torque adds the boolean parameter in the vs.

[ Info: vs = Any[false, 1, 1, 1000, 0.0, 10, [0.12043084138093132 -1.1140463077640161 0.5778026986504824 0.07683217428721806], 1, 1, 0, 1, 1, 2]

In all other cases so far convert(, ) would return a 1.0 or 0.0 prompting no errors. Note that this issue has popped up here because we have array in vs; thus C is promoted to Union{Float64, Int64, Matrix{Float64}} - here^1.

This isn't an issue with MSL@1, as use_support wasn't a parameter.

baggepinnen commented 1 year ago

Ideally, use_support should not be included as a parameter. If I after the creation of the model change use_support before calling solve, will I get the other branch in this code? https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/main/src/Mechanical/Rotational/utils.jl#L142C1-L149C8 If not, I have now introduced a silent bug where I think that I'm running with support, but the model was created without and is thus incorrect. If I indeed do get the other branch, then I'm fine with this, but since the other branch changes the number of equations and variables I'm not sure how that would work

YingboMa commented 1 year ago

This is fixed on the latest MTK and MTKStdlib.