-
-
Notifications
You must be signed in to change notification settings - Fork 209
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Linearization of improper (differentiating / inverse) models #3057
Comments
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 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
What you can do right now is to introduce new input variables 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 = :io
in MRE) works. Likewise, linearization (symbolic and numeric) of the input-output model works (case = :lin_io
in MRE). Also simulation of the zero dynamics works (specifying the outputs;case = :zero
in MRE).HOWEVER, both symbolic and numeric linearization fails for the zero dynamics case (
case = :lin_zero
in 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.
The text was updated successfully, but these errors were encountered: