baggepinnen
commented
1 month ago The generated dynamics function appears to contain unexpanded derivatives
uf.f.f.f_oop
RuntimeGeneratedFunction(#=in ModelingToolkit=#, #=using ModelingToolkit=#, :((ˍ₋arg1, ˍ₋arg2, t)->begin
begin
var"h_1(t)" = (/)(ˍ₋arg2[9], ˍ₋arg2[13])
var"h_1ˍt(t)" = (/)((Differential(t))(ˍ₋arg2[9]), ˍ₋arg2[13])
...The problem is in simplification:
io = ([t4_LZ.y_1, t4_LZ.y_2], [t4_LZ.h_3, t4_LZ.h_4])
ssys, _ = structural_simplify(t4_LZ, io)
julia> full_equations(ssys)
2-element Vector{Equation}:
Differential(t)(h_3(t)) ~ (-a_3*sqrt(2g*h_3(t))) / A_3 + (A_2*((a_2*sqrt((2g*y_2(t)) / k_c)) / A_2 + (-a_4*sqrt(2g*h_4(t))) / A_2 + Differential(t)(y_2(t)) / k_c)*(1 - γ_2)) / (A_3*γ_2)
Differential(t)(h_4(t)) ~ (-a_4*sqrt(2g*h_4(t))) / A_4 + (A_1*((-a_3*sqrt(2g*h_3(t))) / A_1 + Differential(t)(y_1(t)) / k_c + (a_1*sqrt((2g*y_1(t)) / k_c)) / A_1)*(1 - γ_1)) / (A_4*γ_1)there should be no Differential(t) in the equations at this stage.
I'm not surprised that this happens, to invert a strictly proper model you will need some order of derivatives of the inputs. When simulating, this may work out if the input function is differentiable, but when linearizing there is nothing to differentiate, so the compiler has to be augmented with support for this use case, i.e., by
- Erroring with a useful error message.
- Or better, indicate that the user must also provide the derivatives of inputs as another set of inputs
What you can do right now is to introduce new input variables yd_1 and define y_1 as D(y_1) ~ yd_1, I think that this should work
EDIT: It does work
@mtkmodel FourTank begin
# Ref: Johansson, K.H. (2000): "The Quadruple-Tank Process: A Multivariable
# Laboratory Process with an Adjustable Zero",
# IEEE Trans. Ctrl. Tech., Vol. 8, No. 3, May 2000
# 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
# Constant
g = 981, [description = "Acceleration of gravity, cm/s^2"]
# Parameters
# -- tank parameters
A_1 = 28, [description = "Tank 1, cross-section area, cm^2"]
A_2 = 32, [description = "Tank 2, cross-section area, cm^2"]
A_3 = 28, [description = "Tank 3, cross-section area, cm^2"]
A_4 = 32, [description = "Tank 4, cross-section area, cm^2"]
a_1 = 0.071, [description = "Tank 1, cross-section area of outlet hole, cm^2"]
a_2 = 0.057, [description = "Tank 2, cross-section area of outlet hole, cm^2"]
a_3 = 0.071, [description = "Tank 3, cross-section area of outlet hole, cm^2"]
a_4 = 0.057, [description = "Tank 4, cross-section area of outlet hole, cm^2"]
# -- valve parameters, case of LHP zeros/stable zero dynamics
γ_1 = 0.70, [description = "Valve 1 split setting, -"]
γ_2 = 0.60, [description = "Valve 2 split setting, -"]
k_1 = 3.33, [description = "Valve 1 gain, cm3/V.s"]
k_2 = 3.35, [description = "Valve 2 gain, cm3/V.s"]
# -- sensor parameter
k_c = 0.5, [description = "Sensor gain, V/cm"]
end
#
# Model variables, with initial values needed
@variables begin
# Variables
# -- states
h_1(t)=12.4, [description = "Tank 1 level, cm"]
h_2(t)=12.7, [description = "Tank 2 level, cm"]
h_3(t)=1.8, [description = "Tank 3 level, cm"]
h_4(t)=1.4, [description = "Tank 4 level, cm"]
# -- sensor signals
y_1(t), [description = "Output 1: level in Tank 1, cm"]
y_2(t), [description = "Output 2: level in Tank 2, cm"]
dy_1(t) = 0, [description = "Derivative of output 1, cm/s"]
dy_2(t) = 0, [description = "Derivative of output 2, cm/s"]
# -- actuator signals
v_1(t), [description = "Valve 1 signal, V"]
v_2(t), [description = "Valve 2 signal, V"]
end
#
# Model equations
@equations begin
# ODEs
Dt(h_1) ~ -a_1/A_1*sqrt(2g*h_1) + a_3/A_1*sqrt(2g*h_3) + γ_1*k_1/A_1*v_1
Dt(h_2) ~ -a_2/A_2*sqrt(2g*h_2) + a_4/A_2*sqrt(2g*h_4) + γ_2*k_2/A_2*v_2
Dt(h_3) ~ -a_3/A_3*sqrt(2g*h_3) + (1-γ_2)*k_2/A_3*v_2
Dt(h_4) ~ -a_4/A_4*sqrt(2g*h_4) + (1-γ_1)*k_1/A_4*v_1
# Outputs
y_1 ~ k_c*h_1
y_2 ~ k_c*h_2
dy_1 ~ Dt(y_1)
dy_2 ~ Dt(y_2)
# Inputs
if case == :io
v_1 ~ u_v_1(t)
v_2 ~ u_v_2(t)
elseif case == :lin_io
#v_1 ~ u_v_1(t)
#v_2 ~ u_v_2(t)
elseif case == :zero
y_1 ~ u_r_1(t)
y_2 ~ u_r_2(t)
elseif case == :lin_zero
#y_1 ~ u_r_1(t)
#y_2 ~ u_r_2(t)
end
end
end
#
# Instantiating case for linearizing zero dynamics
@named t4_LZ= FourTank(; case=:lin_zero)
t4_LZ= complete(t4_LZ)
io = [t4_LZ.dy_1, t4_LZ.dy_2], [t4_LZ.h_3, t4_LZ.h_4]
#
# Symbolic linearization -- gives error
mats,t4z_ = ModelingToolkit.linearize_symbolic(t4_LZ, io...)
#
# Numeric linearization -- gives error
mats,t4z_ = ModelingToolkit.linearize(t4_LZ, io...;
op=Dict(t4_LZ.y_1 => 6.2, t4_LZ.y_2 => 6.35))
Describe the bug 🐞
Consider the well-known, academic "Four Tank" model. Simulation of the input-output (valve-signal to output levels,
case = :ioin MRE) works. Likewise, linearization (symbolic and numeric) of the input-output model works (case = :lin_ioin MRE). Also simulation of the zero dynamics works (specifying the outputs;case = :zeroin MRE).HOWEVER, both symbolic and numeric linearization fails for the zero dynamics case (
case = :lin_zeroin MRE).Expected behavior
Both symbolic and numeric linearization of the zero dynamics case (
case = :lin_zero) should work. The eigenvalues of the zero dynamics system matrix should equal the transmission zeros ofcase = :lin_io.Minimal Reproducible Example 👇
Here is my case, set up with both symbolic and numeric linearization of the zero dynamics.
Error & Stacktrace ⚠️ Case of symbolic linearization
Case of (numeric) linearization:
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
Add any other context about the problem here.