@@ -27,20 +27,17 @@ rc = 0.25 # Reference concentration
27
27
k0 = 1.05e14
28
28
ϵ = 34.2894
29
29
end
30
-
31
30
@variables begin
32
31
gamma (t), [description = " Reaction speed" ]
33
32
xc (t) = c0, [description = " Concentration" ]
34
33
xT (t) = T0, [description = " Temperature" ]
35
34
xT_c (t), [description = " Cooling temperature" ]
36
35
end
37
-
38
36
@components begin
39
37
T_c = RealInput ()
40
38
c = RealOutput ()
41
39
T = RealOutput ()
42
40
end
43
-
44
41
begin
45
42
τ0 = 60
46
43
wk0 = k0 / c0
@@ -57,20 +54,18 @@ rc = 0.25 # Reference concentration
57
54
gamma ~ xc * wk0 * exp (- wϵ / xT)
58
55
D (xc) ~ - wa11 * xc - wa12 * gamma + wa13
59
56
D (xT) ~ - wa21 * xT + wa22 * gamma + wa23 + wb * xT_c
60
-
61
57
xc ~ c. u
62
58
xT ~ T. u
63
59
xT_c ~ T_c. u
64
60
end
65
61
end
66
-
67
62
begin
68
- Ftf = tf (1 , [(100 ), 1 ])^ 3
63
+ Ftf = tf (1 , [(100 ), 1 ])^ 2
69
64
Fss = ss (Ftf)
70
-
71
65
# Create an MTK-compatible constructor
72
66
function RefFilter (; name)
73
67
sys = ODESystem (Fss; name)
68
+ " Compute initial state that yields y0 as output"
74
69
empty! (ModelingToolkit. get_defaults (sys))
75
70
return sys
76
71
end
82
77
x10 = 0.42
83
78
x20 = 0.01
84
79
u0 = - 0.0224
85
-
86
80
c_start = c0 * (1 - x10) # Initial concentration
87
81
T_start = T0 * (1 + x20) # Initial temperature
88
82
c_high_start = c0 * (1 - 0.72 ) # Reference concentration
104
98
@equations begin
105
99
connect (ref. output, :r , filter. input)
106
100
connect (filter. output, inverse_tank. c)
107
-
108
101
connect (inverse_tank. T_c, ff_gain. input)
109
102
connect (ff_gain. output, :uff , limiter. input)
110
103
connect (limiter. output, add. input1)
111
-
112
104
connect (controller. ctr_output, :u , add. input2)
113
-
114
- # connect(add.output, :u_tot, limiter.input)
115
- # connect(limiter.output, :v, tank.T_c)
116
-
117
105
connect (add. output, :u_tot , tank. T_c)
118
-
119
106
connect (inverse_tank. T, feedback. input1)
120
-
121
107
connect (tank. T, :y , noise_filter. input)
122
-
123
108
connect (noise_filter. output, feedback. input2)
124
109
connect (feedback. output, :e , controller. err_input)
125
110
end
@@ -128,8 +113,10 @@ end;
128
113
ssys = structural_simplify (model)
129
114
cm = complete (model)
130
115
131
- op = Dict (D (cm. inverse_tank. xT) => 1 ,
132
- cm. tank. xc => 0.65 )
116
+ op = Dict (
117
+ cm. filter. y. u => 0.8 * (1 - 0.42 ),
118
+ D (cm. filter. y. u) => 0
119
+ )
133
120
tspan = (0.0 , 1000.0 )
134
121
# https://github.com/SciML/ModelingToolkit.jl/issues/2786
135
122
prob = ODEProblem (ssys, op, tspan)
@@ -151,37 +138,34 @@ p = parameter_values(Sf)
151
138
# https://github.com/SciML/ModelingToolkit.jl/issues/2786
152
139
matrices1 = Sf (x, p, 0 )
153
140
matrices2, _ = Blocks. get_sensitivity (model, :y ; op); # Test that we get the same result when calling the higher-level API
154
- @test matrices1. f_x ≈ matrices2. A[1 : 7 , 1 : 7 ]
141
+ @test matrices1. f_x ≈ matrices2. A[1 : 6 , 1 : 6 ]
155
142
nsys = get_named_sensitivity (model, :y ; op) # Test that we get the same result when calling an even higher-level API
156
143
@test matrices2. A ≈ nsys. A
157
144
158
145
# Test the same thing for comp sensitivities
159
146
160
147
# This should work without providing an operating opint containing a dummy derivative
161
- Sf, simplified_sys = Blocks. get_comp_sensitivity_function (
162
- model, :y ; op, initialization_solver_alg = NewtonRaphson ());
148
+ Sf, simplified_sys = Blocks. get_comp_sensitivity_function (model, :y ; op);
163
149
x = state_values (Sf)
164
150
p = parameter_values (Sf)
165
151
# If this somehow passes, mention it on
166
152
# https://github.com/SciML/ModelingToolkit.jl/issues/2786
167
153
matrices1 = Sf (x, p, 0 )
168
154
# Test that we get the same result when calling the higher-level API
169
- matrices2, _ = Blocks. get_comp_sensitivity (
170
- model, :y ; op, initialization_solver_alg = NewtonRaphson ())
171
- @test matrices1. f_x ≈ matrices2. A[1 : 7 , 1 : 7 ]
155
+ matrices2, _ = Blocks. get_comp_sensitivity (model, :y ; op)
156
+ @test matrices1. f_x ≈ matrices2. A[1 : 6 , 1 : 6 ]
172
157
# Test that we get the same result when calling an even higher-level API
173
- nsys = get_named_comp_sensitivity (
174
- model, :y ; op, initialization_solver_alg = NewtonRaphson ())
158
+ nsys = get_named_comp_sensitivity (model, :y ; op)
175
159
@test matrices2. A ≈ nsys. A
176
160
177
161
@testset " Issue #3319" begin
178
162
op1 = Dict (
179
- D ( cm. inverse_tank . xT) => 1 ,
163
+ cm. filter . y . u => 0.8 * ( 1 - 0.42 ) ,
180
164
cm. tank. xc => 0.65
181
165
)
182
166
183
167
op2 = Dict (
184
- D ( cm. inverse_tank . xT) => 1 ,
168
+ cm. filter . y . u => 0.8 * ( 1 - 0.42 ) ,
185
169
cm. tank. xc => 0.45
186
170
)
187
171
0 commit comments