I’m simulating two nonlinear RC circuits in series (nonlinear resistances, ipp(u) instead of ipp=u/R).
ModelingToolkit code + simulation
# 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.:
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
@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)
This transfer function indicates:
- One time constant at ca. 0.7 seconds
- One time constant at ca. 2.35\cdot 10^6 seconds (almost an integrator)
- One zero at ca. 1.4 seconds
Questions:
- Is there something wrong in my linearization procedure?
- 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??).