MathematicalSystems.jl icon indicating copy to clipboard operation
MathematicalSystems.jl copied to clipboard

IVP macro should return affine system

Open mforets opened this issue 4 years ago • 14 comments

julia> b = [1.0]
1-element Array{Float64,1}:
 1.0

julia> @system(x' = x + b)
LinearControlContinuousSystem{Float64,IdentityMultiple{Float64},IdentityMultiple{Float64}}([1.0], [1.0])

The answer is OK but I would expect an AffineContinuousSystem instead. This way we can use the type information to know that the affine term is a singleton (and not Bu(t) as in the LCCS which is much more general).

mforets avatar Dec 02 '20 13:12 mforets

This at least works

@system(x' = x + b, input=u)
AffineContinuousSystem{Float64,IdentityMultiple{Float64},Array{Float64,1}}([1.0], [1.0])

So it has to interpret b as input, but I am not sure where this happens yet. I'll try to have a look.

ueliwechsler avatar Dec 03 '20 10:12 ueliwechsler

Ah wow a really interesting "feature" of the implementation 😅

julia> @capture(:(x = x + b),  (x_ = f_(x_, u_)))
true
julia> f
:+

So maybe just add a condition here to not allow f==:+ and f==:-? https://github.com/JuliaReach/MathematicalSystems.jl/blob/80abafbc72829b6c82ac3fce06357b273e9337fb/src/macros.jl#L285-L291

ueliwechsler avatar Dec 03 '20 10:12 ueliwechsler

Very neat! Yesterday I thought it would be more involved, but you are right, maybe we just have to add one new condition in line 289, ie. be sure that @system(x' = x + b) does not fallback to the case of an input u.

mforets avatar Dec 03 '20 10:12 mforets

I am surprised that this case is not captured by any test. If it is not that urgent, I can try to make a PR on the weekend.

  • catch f==:+ || f==:- on line 289
  • update the unit tests with @system(x' = x + b) and @system(x' = x - b)

Are there any other corner cases that are not working? Or should be added to the tests?

ueliwechsler avatar Dec 03 '20 10:12 ueliwechsler

Nice! There is no rush. I think adding unit tests as you propose is more than enough.

Are there any other corner cases that are not working?

Do you mean for this issue? (I think the answer is no, at least we haven't found more cases).

On the other hand, there are some open issues related to the macros. I will now create a new label macros so it is easier to find them in the issue tracker.

mforets avatar Dec 03 '20 11:12 mforets

Do you mean for this issue?

Exactly.

On the other hand, there are some open issues related to the macros. I will now create a new label macros so it is easier to find them in the issue tracker.

I am (was) working on these issues on my holidays. There are two branches still open, that I have to make a PR from, that solve some of the issues. I hope to come back to them once I have a bit more time.

ueliwechsler avatar Dec 03 '20 11:12 ueliwechsler

The fix works for the example. But every time I tried to fix something, I figure out that there are more corner cases that are not caught by the logic or ar have not a clean implementation.

https://github.com/JuliaReach/MathematicalSystems.jl/blob/80abafbc72829b6c82ac3fce06357b273e9337fb/src/macros.jl#L411-L414

For example here, we do not consider minus signs which is why the following does not work:

julia> @system(x' = x - [1.0])
ERROR: ArgumentError: the entry () does not match a `MathematicalSystems.jl` structure

Additionally, the expression :(a - b + c) is parsed as true by @capture(rhs, A_ + B__) where A=:(a-b) and B=:c Therefore, it is really hard to fix that with how @system is structured 😐

ueliwechsler avatar Dec 04 '20 01:12 ueliwechsler

The @system (and @ivp) macros have been super useful and we have used them extensively (see for example the ReachabilityModels.jl library). OTOH I have often found that it is not very easy to fix or to generalize for new use cases, such as the new second order system types.

Maybe we should think about the inner workings of the macros without breaking the existing coverage, with the aim of making them easier to fix corner cases and to create new use cases. I'm not sure how, but an idiomatic way may be to have a mechanism that lets you define system's patterns (*) and a dictionary of registered patterns to which to compare the expression given by the user.

(*) for example x' = A*x, x' = A*x + B*u, x' = A*x + B*u + c, or x' = A*x + c. In some cases there may be more than one possible match, so we have to define a preference rule.

mforets avatar Dec 04 '20 11:12 mforets

The idea is that if I want that @system(M*x'' + C*x' + K*x = 0) define a new SecondOrderLinearContinuousSystem, it is sufficient to do something similar to:

register(expr=:(M*x'' + C*x' + K*x = 0), vars=:x, fields=(:M, :C, :K), target=SecondOrderLinearContinuousSystem)
register(expr=:(M*x'' + K*x = 0), vars=:x, fields=(:M, :zeros(...), :K), target=SecondOrderLinearContinuousSystem)
# other cases...

mforets avatar Dec 04 '20 11:12 mforets

I agree with you. It is too much work to extent / adjust the current source code of @system. I think a big issue is, that we wanted to be as generic as possible (e.g. allowing for equations with * and without, let the macro figure out what variables are a state or note, etc...).

I like your proposed idea, it looks quite clean and extendable. However, it is not clear to me, if it is easier to implement this. We might have the same problem in the end. Maybe the flexibility of the macro should be limited. E.g. the equations need to have * and there is only a set of variables that define what is parsed as a state or not. Maybe also with a register mechanism if the user wants to extend it.

I am afraid, I will most likely not have a lot of time to help to refactor the macro in the next months.

But to add to your idea: What about having a check at the beginning of the @system macro, where we check e.g. for a registered second-order system (if this is possible with all the ' since they are parsed rather ugly) and so we could gradually exchange the functionality.

ueliwechsler avatar Dec 04 '20 21:12 ueliwechsler

Concerning your idea and the second order system, I played a bit around and it might be not too difficult to do it, if we can figure out how to parse expressions to the @capture macro mask and put everything into a function (see code below)

I did achieve the following:

julia> @system(M*x'' + C*x' + K*x = 0)
SecondOrderLinearContinuousSystem{Float64,Array{Float64,2},Array{Float64,2},Array{Float64,2}}([1.0 2.0; 2.0 3.0], [1.0 2.0; 2.0 3.0], [1.0 2.0; 2.0 3.0])

with the idea outlined in the post before.

See the source code here: https://github.com/JuliaReach/MathematicalSystems.jl/blob/4bce55c751d90a1222ef245b4ea38443c5c1eee9/src/macros.jl#L886-L889

ueliwechsler avatar Dec 04 '20 23:12 ueliwechsler

(e.g. allowing for equations with * and without,

If this is a problem I'm fine with removing that feature (or keep it but in some other way, like with a global option).

mforets avatar Dec 08 '20 12:12 mforets

Reopening as there are some ideas left in this thread. The second order system construction is nice, I opened a draft PR https://github.com/JuliaReach/MathematicalSystems.jl/pull/225 . I'll play with it.

mforets avatar Dec 08 '20 12:12 mforets

(e.g. allowing for equations with * and without,

If this is a problem I'm fine with removing that feature (or keep it but in some other way, like with a global option).

Even though I like the feature, I think it is not worth to keep it around in a second iteration. Since it introduces "non-mathematical" parsing pattern and therefore does not generalize well.

ueliwechsler avatar Dec 09 '20 17:12 ueliwechsler