# Linearization & ModelingToolkit -- what is wrong?

**URL:** <https://discourse.julialang.org/t/linearization-modelingtoolkit-what-is-wrong/122766>\
**Category:** Modelling & Simulations\
**Created:** [November 18, 2024, 11:36am UTC](https://discourse.julialang.org/t/linearization-modelingtoolkit-what-is-wrong/122766 "2024-11-18T11:36:33Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 18, 2024, 11:36am UTC](https://discourse.julialang.org/t/linearization-modelingtoolkit-what-is-wrong/122766/1 "2024-11-18T11:36:34Z")

</div>

I’m simulating two nonlinear RC circuits in series (nonlinear resistances, `ipp(u)` instead of `ipp=u/R`).

> **ModelingToolkit code + simulation**
>
> ```julia
> # Importing packages
> using ModelingToolkit
> using ModelingToolkit: t_nounits as t, D_nounits as Dt
> using OrdinaryDiffEq
> using ControlSystems, ControlSystemsMTK
> using Plots, LaTeXStrings
> #
> # Nonlinear voltage-> current relationships
> function ipp_1(u)
> R = 8.314
> T = 273.15+60
> F = 9.649e4
> E_a_1 = 18e3
> ipp_s_1 = 0.4607
> nu_1 = 2
> ipp_0_1 = ipp_s_1*exp(-E_a_1/R/T)
> return 2*ipp_0_1*sinh(u*nu_1/2*F/(R*T))
> end
> #
> function ipp_2(u)
> R = 8.314
> T = 273.15+60
> F = 9.649e4
> E_a_2 = 76e3
> ipp_s_2 = 176.1
> nu_2 = 2
> ipp_0_2 = ipp_s_2*exp(-E_a_2/R/T)
> return 2*ipp_0_2*sinh(u*nu_2/2*F/(R*T))
> end
> #
> # Inputs
> @register_symbolic u_ipp(t)
> u_ipp_const(t) = 2
> u_ipp_step(t) = t < 2e-3 ? 2 : 1.9
> #
> # ModelingToolkit model
> @mtkmodel Circuit begin
> #
> # System structure parameters
> @structural_parameters begin
> # :io = standard input-output, :zero = zero dynamics
> # :lin_io = linearizing input-output, :lin_zero = linearizing zero dynamics
> case = :io
> end
> #
> # Model parameters
> @parameters begin
> # Parameters
> A = 160, [description = "Area, cm2"]
> Cpp = 350e-4, [description = "Capacitance per unit area, F/cm^2"]
> C = A*Cpp, [description = "Capacitance, F"]
> delta_m = 100e-6*100, [description = "Thickness, cm"]
>        
> end
> #
> # Model variables, with initial values needed
> @variables begin
> # Variables
> i(t), [description = "Current, A"]
> i_1(t)
> i_2(t)
> i_r_1(t), [description = "Resistor 1 current, A"]
> i_c_1(t), [description = "Capacitor 1 current, A"]
> u_1(t)=0, [description = "Voltage 1, V"]
> i_r_2(t), [description = "Resistor 2 urrent, A"]
> i_c_2(t), [description = "Capacitor 2 current, A"]
> u_2(t)=0, [description = "Voltage 2, V"]
> u(t), [description = "Total voltage, V"]
> #
> # Control inputs
> T(t), [description = "System temperature, K"]
> ipp(t), [description = "Current density, A/cm2"] 
> end
> #
> # Model equations
> @equations begin
> # Balance differential equations
> i ~ ipp*A
> i_r_1 ~ A*ipp_1(u_1)
> i_1 ~ i_r_1 + i_c_1
> Dt(u_1) ~ i_c_1/C
> i_r_2 ~ A*ipp_2(u_2)
> i_2 ~ i_r_2 + i_c_2
> Dt(u_2) ~ i_c_2/C
> u ~ u_1 + u_2
> i ~ i_1 
> i ~ i_2 
> # Optionally leave out input when model is used for linearization
> if case == :io
> ipp ~ u_ipp(t)
> end
> end
> end
> #
> # Instantiating and simulating system
> @mtkbuild cir = Circuit()
> tspan = (0,1)
> # -- simulate with constant input until steady state
> u_ipp(t) = u_ipp_const(t)
> prob = ODEProblem(cir, [], tspan)
> sol = solve(prob, Rodas5(); reltol=1e-5,tstops=[20*60])
> u_ss = sol(tspan[2])
> # -- restart new simulation at steady state
> tspan = (0,1e-2)
> u_ipp(t) = u_ipp_step(t)
> prob = ODEProblem(cir, u_ss, tspan)
> sol = solve(prob, Rodas5(); reltol=1e-6,tstops=[0,2e-3])
> plot(sol, idxs=(t*1e3,cir.u), lw=2.5, label=L"u", xlims=(0,10), xlabel=L"$t$ [ms]")
> 
> ```

Here is a simple step response of the total voltage across the two circuits, using ModelingToolkit, etc.:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/a/7ac7d41ec094dc9752005c058b227ddfe51f5d85.png)

Using the “63% rule”, the time constant should be in the order of 0.5 _milliseconds_, or so [NOTE: I’ve multiplied the time axis by `1e3` in the plot].

However, if I linearize the model using ModelingToolkit/ControlSystemsMTK, I get the following transfer function:

> **ModelingToolkit linearization**
>
> ```julia
> @named cir_LIO = Circuit(; case = :lin_io)
> cir_LIO = complete(cir_LIO)
> sys = named_ss(cir_LIO, [cir_LIO.ipp], [cir_LIO.u], op=Dict(cir_LIO.ipp=>2.0))
> zpk(sys)
> 
> ```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/8/08a95941b7b9fcb08d3bc6285d804e5364f2b6c3.png)  
This transfer function indicates:

1. One time constant at ca. 0.7 _seconds_
2. One time constant at ca. 2.35\cdot 10^6 seconds (almost an integrator)
3. One zero at ca. 1.4 seconds

**Questions** :

1. Is there something wrong in my linearization procedure?
2. Is there a numeric problem with this system? (The system is very low order, so even if it is nonlinear, it should be do-able??).

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 18, 2024, 12:33pm UTC](https://discourse.julialang.org/t/linearization-modelingtoolkit-what-is-wrong/122766/2 "2024-11-18T12:33:56Z")

</div>

Your operating point does not include non-zero values for all variables you may be interested in, with

```julia
op = Dict(unknowns(cir) .=> u_ss)
op = merge(op, Dict(cir_LIO.ipp=>2.0))
sys = named_ss(cir_LIO, [cir_LIO.ipp], [cir_LIO.u]; op)

```

I get

```julia
julia> zpk(sys)
TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
                               1.0s + 1990.6480012656318
57.14285714285714------------------------------------------------------
                 (1.0s + 1990.6482406641012)(1.0s + 1990.6477618671631)

Continuous-time transfer function model

julia> dampreport(sys)
| Pole | Damping | Frequency | Frequency | Time Constant |
| | Ratio | (rad/sec) | (Hz) | (sec) |
+--------------------+---------------+---------------+---------------+---------------+
| -1.99e+03 | 1 | 1.99e+03 | 317 | 0.000502 |
| -1.99e+03 | 1 | 1.99e+03 | 317 | 0.000502 |

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 18, 2024, 1:15pm UTC](https://discourse.julialang.org/t/linearization-modelingtoolkit-what-is-wrong/122766/3 "2024-11-18T13:15:18Z")

</div>

> [@baggepinnen](#):
>
> `op = Dict(unknowns(cir) .=> u_ss)`

Ah. I mistakenly thought that I had used `u_ss`, but I see that I only did that for the simulation without linearization – and that the linearized model did not know about this.

Thanks!