permute inputs in inputs_to_parameters!

to correspond to user-provided input ordering. Otherwise subsequent linear algebra operations will make no sense.

Author: baggepinnen

src/inputoutput.jl

@@ -218,7 +218,8 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     f, dvs, ps
 end
 
-function inputs_to_parameters!(state::TransformationState, check_bound = true)
+function inputs_to_parameters!(state::TransformationState, io)
+    check_bound = io === nothing
     @unpack structure, fullvars, sys = state
     @unpack var_to_diff, graph, solvable_graph = structure
     @assert solvable_graph === nothing
@@ -274,6 +275,15 @@ function inputs_to_parameters!(state::TransformationState, check_bound = true)
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))
     ps = parameters(sys)
+
+    if io !== nothing
+        # Change order of new parameters to correspond to user-provided order in argument `inputs`
+        param_permutation = map(io.inputs) do inp
+            findfirst(isequal(inp), new_parameters)
+        end
+        new_parameters = new_parameters[param_permutation]
+    end
+
     @set! sys.ps = [ps; new_parameters]
 
     @set! state.sys = sys

src/systems/abstractsystem.jl

@@ -965,7 +965,7 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     state = TearingState(sys)
     has_io = io !== nothing
     has_io && markio!(state, io...)
-    state, input_idxs = inputs_to_parameters!(state, !has_io)
+    state, input_idxs = inputs_to_parameters!(state, io)
     sys = alias_elimination!(state)
     # TODO: avoid construct `TearingState` again.
     state = TearingState(sys)
@@ -981,7 +981,7 @@ end
 
 function io_preprocessing(sys::AbstractSystem, inputs,
                           outputs; simplify = false, kwargs...)
-    sys, input_idxs = structural_simplify(sys, (inputs, outputs); simplify, kwargs...)
+    sys, input_idxs = structural_simplify(sys, (; inputs, outputs); simplify, kwargs...)
 
     eqs = equations(sys)
     alg_start_idx = findfirst(!isdiffeq, eqs)
@@ -1215,7 +1215,7 @@ function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives =
         A = [f_x f_z
              gzgx*f_x gzgx*f_z]
         B = [f_u
-             gzgx*f_u]
+             gzgx * f_u] # The cited paper has zeros in the bottom block, see derivation in https://github.com/SciML/ModelingToolkit.jl/pull/1691 for the correct formula
 
         C = [h_x h_z]
         Bs = -(gz \ g_u) # This equation differ from the cited paper, the paper is likely wrong since their equaiton leads to a dimension mismatch.