refactor: directly solve initialization problem in `linearization_function` #2762

AayushSabharwal · 2024-06-04T06:29:07Z

Checklist

Appropriate tests were added
Any code changes were done in a way that does not break public API
All documentation related to code changes were updated
The new code follows the
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.
Any new documentation only uses public API

Additional context

Add any other context about the problem here.

baggepinnen · 2024-06-04T09:26:58Z

src/systems/abstractsystem.jl

@@ -1749,6 +1749,7 @@ The `simplified_sys` has undergone [`structural_simplify`](@ref) and had any occ
  - `outputs`: A vector of variables that indicate the outputs of the linearized input-output model.
  - `simplify`: Apply simplification in tearing.
  - `initialize`: If true, a check is performed to ensure that the operating point is consistent (satisfies algebraic equations). If the op is not consistent, initialization is performed.
+  - `initialization_solver_alg`: A NonlinearSolve algorithm to use for solving the initialization problem.


Suggested change

- `initialization_solver_alg`: A NonlinearSolve algorithm to use for solving the initialization problem.

- `initialization_solver_alg`: A NonlinearSolve algorithm to use for solving for a feasible set of state and algebraic variables that satisfies the specified operating point.

"initialization" has no well-defined meaning when linearizing, hence the slightly different language.

Fixed in 406dc9a

ChrisRackauckas · 2024-06-05T23:15:34Z

Why is this just a draft?

AayushSabharwal · 2024-06-06T05:13:47Z

I wanted to make sure it's rebased before merging

baggepinnen · 2024-06-10T10:20:27Z

On this PR branch, I get a funky error for this example

using ModelingToolkit, ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks
@parameters t;
D = Differential(t);
rc = 0.25 # Reference concentration 

@mtkmodel MixingTank begin
    @parameters begin
        c0  = 0.8,   [description = "Nominal concentration"]
        T0  = 308.5, [description = "Nominal temperature"]
        a1  = 0.2674
        a21 = 1.815
        a22 = 0.4682
        b   = 1.5476
        k0  = 1.05e14
        ϵ   = 34.2894
    end

    @variables begin
        gamma(t),     [description = "Reaction speed"]
        xc(t) = c0,   [description = "Concentration"]
        xT(t) = T0,   [description = "Temperature"]
        xT_c(t) = T0, [description = "Cooling temperature"]
    end

    @components begin
        T_c = RealInput()
        c = RealOutput()
        T = RealOutput()
    end

    begin
        τ0   = 60
        wk0  = k0/c0
        wϵ   = ϵ*T0
        wa11 = a1/τ0
        wa12 = c0/τ0
        wa13 = c0*a1/τ0
        wa21 = a21/τ0
        wa22 = a22*T0/τ0
        wa23 = T0*(a21 - b)/τ0
        wb   = b/τ0
    end
    @equations begin
        gamma ~ xc*wk0*exp( -wϵ/xT)
        D(xc) ~ -wa11*xc - wa12*gamma + wa13
        D(xT) ~ -wa21*xT + wa22*gamma + wa23 + wb*xT_c

        xc ~ c.u
        xT ~ T.u
        xT_c ~ T_c.u
    end

end


@mtkmodel InverseControlledTank begin
	begin
		c0 = 0.8    #  "Nominal concentration
		T0 = 308.5 	#  "Nominal temperature
		x10 = 0.42 	
		x20 = 0.01 
		u0 = -0.0224 

		c_start = c0*(1-x10) 		# Initial concentration
		T_start = T0*(1+x20) 		# Initial temperature
		c_high_start = c0*(1-0.72) 	# Reference concentration
		T_c_start = T0*(1+u0) 		# Initial cooling temperature
	end
	@components begin
		ref 			= Constant(k=0.25) # Concentration reference
		ff_gain 		= Gain(k=1) # To allow turning ff off
		controller 		= PI(gainPI.k=10, T=500)
		tank 			= MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
		inverse_tank 	= MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
		feedback 		= Feedback()
		add 			= Add()
		filter 			= SecondOrder(k=1, w=1) # Initialize filter states to the initial concentration
		noise_filter 	= FirstOrder(k=1, T=1, x=T_start)
		limiter 		= Limiter(y_max=370, y_min=250) # Saturate the control input 
	end
	@equations begin
		connect(ref.output, :r, filter.input)
		connect(filter.output, inverse_tank.c)
		connect(inverse_tank.T_c, ff_gain.input)
		connect(ff_gain.output, :uff, limiter.input)
		connect(limiter.output, add.input1)
		connect(controller.ctr_output, :u, add.input2)
		connect(add.output, :u_tot, tank.T_c)
		connect(inverse_tank.T, feedback.input1)
		connect(tank.T, :y, noise_filter.input)
		connect(noise_filter.output, feedback.input2)
		connect(feedback.output, :e, controller.err_input)
	end
end;
@named model = InverseControlledTank()
get_sensitivity(model, :y)

julia> get_sensitivity(model, :y)
ERROR: MethodError: no method matching (::ModelingToolkit.var"#fn#209"{…})(::Vector{…}, ::ModelingToolkit.MTKParameters{…})

Closest candidates are:
  (::ModelingToolkit.var"#fn#209")(::Any, ::ModelingToolkit.MTKParameters, ::Any)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/abstractsystem.jl:500
  (::ModelingToolkit.var"#fn#209")(::Any, ::Any, ::Any)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/abstractsystem.jl:499

Stacktrace:
 [1] (::SymbolicIndexingInterface.TimeIndependentObservedFunction{…})(::SymbolicIndexingInterface.NotTimeseries, prob::SciMLBase.NonlinearSolution{…})
   @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/A1VUA/src/state_indexing.jl:96
 [2] (::SymbolicIndexingInterface.TimeIndependentObservedFunction{…})(prob::SciMLBase.NonlinearSolution{…})

AayushSabharwal · 2024-06-10T10:24:27Z

Fixed

baggepinnen · 2024-06-10T12:02:20Z

Thanks for quick fix 🎉 With this slight modification to the example, I get another error

using ModelingToolkit, ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks
using ControlSystemsMTK
@parameters t;
D = Differential(t);
rc = 0.25 # Reference concentration 

@mtkmodel MixingTank begin
    @parameters begin
        c0  = 0.8,   [description = "Nominal concentration"]
        T0  = 308.5, [description = "Nominal temperature"]
        a1  = 0.2674
        a21 = 1.815
        a22 = 0.4682
        b   = 1.5476
        k0  = 1.05e14
        ϵ   = 34.2894
    end

    @variables begin
        gamma(t),     [description = "Reaction speed"]
        xc(t) = c0,   [description = "Concentration"]
        xT(t) = T0,   [description = "Temperature"]
        xT_c(t) = T0, [description = "Cooling temperature"]
    end

    @components begin
        T_c = RealInput()
        c = RealOutput()
        T = RealOutput()
    end

    begin
        τ0   = 60
        wk0  = k0/c0
        wϵ   = ϵ*T0
        wa11 = a1/τ0
        wa12 = c0/τ0
        wa13 = c0*a1/τ0
        wa21 = a21/τ0
        wa22 = a22*T0/τ0
        wa23 = T0*(a21 - b)/τ0
        wb   = b/τ0
    end
    @equations begin
        gamma ~ xc*wk0*exp( -wϵ/xT)
        D(xc) ~ -wa11*xc - wa12*gamma + wa13
        D(xT) ~ -wa21*xT + wa22*gamma + wa23 + wb*xT_c

        xc ~ c.u
        xT ~ T.u
        xT_c ~ T_c.u
    end

end

begin
	Ftf = tf(1, [(100), 1])^3
	Fss = ss(Ftf)

	"Compute initial state that yields y0 as output"
	function init_filter(y0)
	    (; A,B,C,D) = Fss
	    Fx0 = -A\B*y0
	    @assert C*Fx0 ≈ [y0] "C*Fx0*y0 ≈ y0 failed, got $(C*Fx0*y0) ≈ $(y0)]" 
	    Fx0
	end

	# Create an MTK-compatible constructor 
	RefFilter(; y0, name) = ODESystem(Fss; name, x0=init_filter(y0))
end
@mtkmodel InverseControlledTank begin
	begin
		c0 = 0.8    #  "Nominal concentration
		T0 = 308.5 	#  "Nominal temperature
		x10 = 0.42 	
		x20 = 0.01 
		u0 = -0.0224 

		c_start = c0*(1-x10) 		# Initial concentration
		T_start = T0*(1+x20) 		# Initial temperature
		c_high_start = c0*(1-0.72) 	# Reference concentration
		T_c_start = T0*(1+u0) 		# Initial cooling temperature
	end
	@components begin
		ref 			= Constant(k=0.25) # Concentration reference
		ff_gain 		= Gain(k=1) # To allow turning ff off
		controller 		= PI(gainPI.k=10, T=500)
		tank 			= MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
		inverse_tank 	= MixingTank(xc=c_start, xT = T_start, c0=c0, T0=T0)
		feedback 		= Feedback()
		add 			= Add()
		filter 			= RefFilter(y0=c_start) # Initialize filter states to the initial concentration
		noise_filter 	= FirstOrder(k=1, T=1, x=T_start)
		limiter 		= Limiter(y_max=370, y_min=250) # Saturate the control input 
	end
	@equations begin
		connect(ref.output, :r, filter.input)
		connect(filter.output, inverse_tank.c)
		connect(inverse_tank.T_c, ff_gain.input)
		connect(ff_gain.output, :uff, limiter.input)
		connect(limiter.output, add.input1)
		connect(controller.ctr_output, :u, add.input2)
		connect(add.output, :u_tot, tank.T_c)
		connect(inverse_tank.T, feedback.input1)
		connect(tank.T, :y, noise_filter.input)
		connect(noise_filter.output, feedback.input2)
		connect(feedback.output, :e, controller.err_input)
	end
end;
@named model = InverseControlledTank()
get_sensitivity(model, :y)

ERROR: UndefVarError: `filter₊x(t)[1]ˍtt` not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/0opve/src/code.jl:375 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [6] f
    @ ~/.julia/packages/ModelingToolkit/3Suzf/src/systems/nonlinear/nonlinearsystem.jl:292 [inlined]
  [7] NonlinearFunction
    @ ~/.julia/packages/SciMLBase/roUH9/src/scimlfunctions.jl:2297 [inlined]
  [8] evaluate_f(prob::NonlinearLeastSquaresProblem{…}, u::Vector{…})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/m8oNC/src/internal/helpers.jl:7
  [9] __init(::NonlinearLeastSquaresProblem{…}, ::NonlinearSolve.GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/m8oNC/src/core/generalized_first_order.jl:164
 [10] __init
    @ ~/.julia/packages/NonlinearSolve/m8oNC/src/core/generalized_first_order.jl:154 [inlined]
 [11] __solve(::NonlinearLeastSquaresProblem{…}, ::NonlinearSolve.GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, kwargs::@Kwargs{})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/m8oNC/src/core/generic.jl:3
 [12] __solve
    @ ~/.julia/packages/NonlinearSolve/m8oNC/src/core/generic.jl:1 [inlined]

AayushSabharwal · 2024-06-10T12:16:20Z

That's the structural_simplify bug

baggepinnen · 2024-06-10T12:17:03Z

Is there an issue associated with this bug?

AayushSabharwal · 2024-06-10T12:18:47Z

There isn't. I can open one.

…ction`

default to `TrustRegion`

baggepinnen reviewed Jun 4, 2024

View reviewed changes

AayushSabharwal marked this pull request as draft June 4, 2024 09:32

AayushSabharwal force-pushed the as/faster-init branch from 8ff1743 to 48c034e Compare June 6, 2024 05:14

AayushSabharwal marked this pull request as ready for review June 6, 2024 05:14

AayushSabharwal force-pushed the as/faster-init branch from 48c034e to b453151 Compare June 10, 2024 10:22

AayushSabharwal mentioned this pull request Jun 10, 2024

providing updated op to linearization_function has no effect #2784

Closed

AayushSabharwal added 4 commits June 12, 2024 11:01

refactor: directly solve initialization problem in `linearization_fun…

b60976d

…ction`

feat: allow specifying initialization solver algorithm

dc395f6

default to `TrustRegion`

fix: respect operating point in linearize

7359ddb

test: mark inversemodel tests as broken

568cd44

AayushSabharwal force-pushed the as/faster-init branch from 5f1893f to 568cd44 Compare June 12, 2024 05:32

ChrisRackauckas merged commit 8295ed1 into SciML:master Jun 12, 2024
20 of 24 checks passed

AayushSabharwal deleted the as/faster-init branch June 13, 2024 05:51

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

refactor: directly solve initialization problem in `linearization_function` #2762

refactor: directly solve initialization problem in `linearization_function` #2762

Uh oh!

AayushSabharwal commented Jun 4, 2024 •

edited

Loading

Uh oh!

baggepinnen Jun 4, 2024

Uh oh!

AayushSabharwal Jun 4, 2024

Uh oh!

ChrisRackauckas commented Jun 5, 2024

Uh oh!

AayushSabharwal commented Jun 6, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024 •

edited

Loading

Uh oh!

Uh oh!

Uh oh!

	- `initialization_solver_alg`: A NonlinearSolve algorithm to use for solving the initialization problem.
	- `initialization_solver_alg`: A NonlinearSolve algorithm to use for solving for a feasible set of state and algebraic variables that satisfies the specified operating point.

Uh oh!

refactor: directly solve initialization problem in linearization_function #2762

refactor: directly solve initialization problem in linearization_function #2762

Uh oh!

Conversation

AayushSabharwal commented Jun 4, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Checklist

Additional context

Uh oh!

baggepinnen Jun 4, 2024

Choose a reason for hiding this comment

Uh oh!

AayushSabharwal Jun 4, 2024

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Jun 5, 2024

Uh oh!

AayushSabharwal commented Jun 6, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024

Uh oh!

baggepinnen commented Jun 10, 2024

Uh oh!

AayushSabharwal commented Jun 10, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

refactor: directly solve initialization problem in `linearization_function` #2762

refactor: directly solve initialization problem in `linearization_function` #2762

AayushSabharwal commented Jun 4, 2024 •

edited

Loading

AayushSabharwal commented Jun 10, 2024 •

edited

Loading