add named_ss method

Author: baggepinnen

src/ode_system.jl

@@ -139,6 +139,47 @@ function RobustAndOptimalControl.named_ss(
     nsys
 end
 
+function RobustAndOptimalControl.named_ss(
+    sys::ModelingToolkit.AbstractSystem, linfun::ModelingToolkit.LinearizationFunction, outputs;
+    descriptor = true,
+    simple_infeigs = true,
+    balance = descriptor && !simple_infeigs, # balance only if descriptor is true and simple_infeigs is false
+    big = false,
+    kwargs...,
+)
+    ssys = sys
+    matrices, xpt = ModelingToolkit.linearize(sys, linfun; kwargs...)
+    inputs = linfun.inputs
+    nu = length(inputs)
+    unames = symstr.(inputs)
+    if nu > 0 && size(matrices.B, 2) == 2nu
+        # This indicates that input derivatives are present
+        duinds = findall(any(!iszero, eachcol(matrices.B[:, nu+1:end]))) .+ nu
+        u2du = (1:nu) .=> duinds # This maps inputs to their derivatives
+        lsys = causal_simplification(matrices, u2du; descriptor, simple_infeigs, big, balance, verbose=false)
+    else
+        lsys = ss(matrices...)
+    end
+    pind = [ModelingToolkit.parameter_index(ssys, i) for i in ModelingToolkit.inputs(ssys)]
+    x0 = xpt.x
+    u0 = [xpt.p[pi] for pi in pind] 
+    xu = (; x = x0, u = u0)
+    extra = Dict(:operating_point => xu)
+    # If simple_infeigs=false, the system might have been reduced and the state names might not match the original system.
+    x_names = get_x_names(lsys, ssys; descriptor, simple_infeigs, balance)
+    nsys = named_ss(
+        lsys;
+        x = x_names,
+        u = unames,
+        y = symstr.(outputs),
+        name = string(Base.nameof(sys)),
+        extra,
+    )
+    RobustAndOptimalControl.set_extra!(nsys, :ssys, ssys)
+    nsys
+
+end
+
 function get_x_names(lsys, sys; descriptor, simple_infeigs, balance)
     generic = if descriptor
         !simple_infeigs || balance
@@ -161,7 +202,7 @@ If `descriptor = true`, the function `DescriptorSystems.dss2ss` is used. In this
 
 The argument `big = true` performs computations in `BigFloat` precision, useful for poorly scaled systems. This may require the user to install and load `GenericLinearAlgebra` (if you get error `no method matching svd!(::Matrix{BigFloat})`).
 """
-function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; balance=false, descriptor=true, simple_infeigs = true, big = false)
+function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; balance=false, descriptor=true, simple_infeigs = true, big = false, verbose = true)
     T = big ? BigFloat : Float64
     b1 = big ? Base.big(1.0) : 1.0
     fm(x) = convert(Matrix{T}, x)
@@ -190,7 +231,7 @@ function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; balance=fa
             dsys, T1, T2 = RobustAndOptimalControl.DescriptorSystems.gprescale(dsys)
         else
             bq = RobustAndOptimalControl.DescriptorSystems.gbalqual(dsys)
-            bq > 10000 && @warn("The numerical balancing of the system is poor (gbalqual = $bq), consider using `balance=true` to balance the system before conversion to StateSpace to improve accuracy of the result.")
+            verbose && bq > 10000 && @warn("The numerical balancing of the system is poor (gbalqual = $bq), consider using `balance=true` to balance the system before conversion to StateSpace to improve accuracy of the result.")
         end
 
         # NOTE: the conversion implemented in ss(dss) uses gss2ss due to it's initial call to gir to produce a reduced order model and then an SVD-based alg to improve numerics. Should we use this by default?