source: branches/dev-5819/Python/src/tests_jmodelica/test_delay.py @ 13461

Last change on this file since 13461 was 13461, checked in by randersson, 3 months ago

#5819 Performed 2to3 on the entire src/Python directory in order to make the code compatible. I've not reviewed all changes while doing this commit but will do while starting to run the tests. I'm not sure everything is converted correctly, for example the change in test_casadi_collocation row 1610.

File size: 15.1 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4# Copyright (C) 2014 Modelon AB
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation, version 3 of the License.
9#
10# This program is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18"""Tests for delay simulation."""
19
20import numpy as N
21
22from tests_jmodelica.general.base_simul import *
23from tests_jmodelica import testattr, get_files_path
24from nose.tools import nottest
25
26
27compiler_options={'compliance_as_warning':True}
28
29path_to_mos = os.path.join(get_files_path(), 'Modelica')
30path_to_delay_mo = os.path.join(path_to_mos, 'TestDelay.mo')
31path_to_fixed_delay_mo = os.path.join(path_to_mos, 'TestFixedDelay.mo')
32path_to_spatialdist_mo = os.path.join(path_to_mos, 'TestSpatialDistribution.mo')
33
34def compile_and_load(class_name, filename=path_to_delay_mo):
35    fmu_name = compile_fmu(class_name, filename, compiler_options=compiler_options)
36    return load_fmu(fmu_name) 
37
38def simulate(fmu, final_time, maxh = None):
39    opts = fmu.simulate_options()     
40    opts['solver'] = 'CVode'           
41    if maxh is not None:
42        opts['CVode_options']['maxh'] = maxh
43
44    res = fmu.simulate(final_time = final_time, options = opts)
45    return res
46
47def compile_and_simulate(class_name, final_time, filename=path_to_delay_mo, maxh = None):
48    fmu = compile_and_load(class_name, filename=filename)
49    return simulate(fmu, final_time, maxh)
50
51def assert_close(x, y, abstol):
52    assert len(x) == len(y)
53    d = N.abs(x-y).max()
54    assert d <= abstol, ("Signals differ by " + str(d) + " which is more than abstol = " + str(abstol) )
55
56def switch_signal(t, yp, ym, s, s_eps = 1e-6):
57    """
58    Return a vector y pointwise choosen from yp an yn trying tp pick yp (ym)
59    when s is positive (negative), switching once when t[k]==t[k+1] and s is small
60    """
61    yp, ym, s = yp+0*t, ym+0*t, s+0*t # make sure they have the same shape as t
62    switches = (N.diff(t) == 0) & (N.abs(s[:-1]) <= s_eps)
63    negative = N.bitwise_xor.accumulate(N.hstack([s[0] < 0, switches]))
64    y = N.array(yp)
65    y[negative] = ym[negative]
66    return y
67
68
69def sol_repeating_events(t, d=1):
70    x_expected = N.mod(t, d)
71    (inds,) = N.nonzero(N.diff(t) == 0)
72    x_expected[inds] = d
73    x_expected[inds+1] = 0
74    return x_expected
75
76
77class TestDelay:
78    @testattr(stddist_full = True)
79    def test_next_event(self):
80        # Test bug where delay in inactive if branch does not compute next time event correctly. #5692.
81        res = compile_and_simulate('TestNextEvent', 2)
82
83
84class TestFixedDelay:
85    def zero_delay_ok(self):
86        return True
87
88    def get_class_postfix(self):
89        return ""
90
91    def compile_and_load(self, class_name):
92        return compile_and_load(class_name+self.get_class_postfix(), filename=path_to_fixed_delay_mo)
93
94    def compile_and_simulate(self, class_name, *args, **kwargs):
95        return compile_and_simulate(class_name+self.get_class_postfix(), *args, filename=path_to_fixed_delay_mo, **kwargs)
96
97    @testattr(stddist_full = True)
98    def test_delay_time(self):
99        res = compile_and_simulate('TestDelayTime', final_time = 5, maxh = 0.5)
100        t, x = res['time'], res['x']
101        x_expected = N.maximum(0, t-1)
102        assert_close(x, x_expected, 1e-8)
103
104    @testattr(stddist_full = True)
105    def test_delay_quadratic(self):
106        fmu = compile_and_load('TestDelayQuadratic')
107        res = simulate(fmu, final_time = 5, maxh = 0.25)
108        t, x = res['time'], res['x']
109        x_expected = N.maximum(0, t-1)**2 + 1
110
111        assert_close(x, x_expected, 1e-4)
112
113        fmu.reset()
114        res = simulate(fmu, final_time = 5, maxh = 1/2.5)
115        t, x = res['time'], res['x']
116        x_expected = N.maximum(0, t-1)**2 + 1
117
118        assert_close(x, x_expected, 0.1)
119
120    @testattr(stddist_full = True)
121    def test_integrate_delayed_time(self):
122        res = self.compile_and_simulate('TestIntegrateDelayedTime', final_time = 5)
123        t, x = res['time'], res['x']
124        x_expected = N.maximum(0, t-1)**2/2
125        assert_close(x, x_expected, 1e-4)
126
127    @testattr(stddist_full = True)
128    def test_integrate_delayed_quadratic(self):
129        res = self.compile_and_simulate('TestIntegrateDelayedQuadratic', final_time = 5)
130        t, x = res['time'], res['x']
131        x_expected = N.maximum(0,t-1)**3/3
132        assert_close(x, x_expected, 0.2)
133
134    @testattr(stddist_full = True)
135    def test_sinusoid(self):
136        res = self.compile_and_simulate('TestSinusoid', final_time = 20)
137        t, x = res['time'], res['x']
138        x_expected = 1.03*N.cos((t+0.35)*N.pi/2)
139        inds = (t >= 0.5)
140        assert_close(x[inds], x_expected[inds], 0.1)
141
142    @testattr(stddist_full = True)
143    def test_sinusoid_noevent(self):
144        # Check that we get the same answer as above with noEvent
145        res = self.compile_and_simulate('TestSinusoidNoEvent', final_time = 20)
146        t, x = res['time'], res['x']
147        x_expected = 1.03*N.cos((t+0.35)*N.pi/2)
148        inds = (t >= 0.5)
149        assert_close(x[inds], x_expected[inds], 0.1)
150
151    @testattr(stddist_full = True)
152    def test_short_delay(self):
153        fmu = self.compile_and_load('TestShortDelay')
154        fmu.set('d', 1e-3)
155        res = simulate(fmu, final_time = 5)
156        t, x = res['time'], res['x']
157        x_expected = N.exp(-t)
158
159        assert_close(x, x_expected, 1e-3)
160        assert len(t) < 200
161
162        if self.zero_delay_ok():
163            fmu.reset()
164            fmu.set('d', 0)
165            res = simulate(fmu, final_time = 5)
166            t, x = res['time'], res['x']
167            x_expected = N.exp(-t)
168
169            assert_close(x, x_expected, 1e-3)
170            assert len(t) < 200
171
172    @testattr(stddist_full = True)
173    def test_commute(self):
174        res = self.compile_and_simulate('TestCommute', final_time = 10, maxh = 1/5.5)
175        t, x_delay, delay_x = res['time'], res['x_delay'], res['delay_x']
176        x_expected = N.cos(N.maximum(0,t-1))
177        assert_close(x_delay, x_expected, 1e-2)
178        assert_close(delay_x, x_expected, 1e-2)
179
180    @testattr(stddist_full = True)
181    def test_repeating_events(self):
182        fmu = self.compile_and_load('TestRepeatingEvents')
183
184        for (k, d) in enumerate([1.0, N.nextafter(1.0, 0), N.nextafter(1.0, 2)]):
185            if k > 0: fmu.reset()
186
187            fmu.set('d', d)
188            res = simulate(fmu, final_time = 5.5, maxh = 1/2.5) 
189            t, x = res['time'], res['x']
190            (inds,) = N.nonzero(N.diff(t) == 0)
191            x_expected = sol_repeating_events(t, d)
192
193            assert_close(t[inds], [1,2,3,4,5], 1e-7)
194            assert_close(x, x_expected, 1e-7)
195
196    #@testattr(stddist_full = True)
197    @nottest # why is it not working?
198    def test_repeat_noevent(self):
199        res = self.compile_and_simulate('TestRepeatNoEvent', final_time = 5.5, maxh = 1/10.5)
200        t, x = res['time'], res['x']
201        (inds,) = N.nonzero(N.diff(t) == 0)
202        x_expected = sol_repeating_events(t)
203
204        assert len(inds) == 1
205        assert N.abs(t[inds]-[1]) <= 1e-8
206        assert sum(N.abs(x-x_expected) > 1e-8) <= 6
207
208class TestFixedDelaySpatialDist(TestFixedDelay):
209    def zero_delay_ok(self):
210        return False
211
212    def get_class_postfix(self):
213        return "(redeclare block FD=FixedDelaySD)"
214
215    # Disable these until we can support noEvent with spatialDistribution and two outputs
216    def test_sinusoid_noevent(self):
217        pass
218    def test_repeat_noevent(self):
219        pass
220
221class TestFixedDelaySpatialDistRev(TestFixedDelay):
222    def zero_delay_ok(self):
223        return False
224
225    def get_class_postfix(self):
226        return "(redeclare block FD=FixedDelaySDReverse)"
227
228@testattr(stddist_full = True)
229def test_repeat_noevent():
230    res = compile_and_simulate('TestRepeatNoEvent', final_time = 5.5, maxh = 1/10.5)
231    t, x = res['time'], res['x']
232    (inds,) = N.nonzero(N.diff(t) == 0)
233    x_expected = sol_repeating_events(t)
234
235    assert len(inds) == 1
236    assert N.abs(t[inds]-[1]) <= 1e-8
237    assert sum(N.abs(x-x_expected) > 1e-8) <= 6
238
239
240@testattr(stddist_full = True)
241def test_variably_delayed_time():
242    res = compile_and_simulate('TestVariablyDelayedTime', final_time = 5, maxh = 0.1)
243    t, x = res['time'], res['x']
244    x_expected = N.maximum(0,t-(N.sin(5*t)*0.5+0.5))
245    assert_close(x, x_expected, 1e-7)
246
247@testattr(stddist_full = True)
248def test_state_dependent_delay_time():
249    res = compile_and_simulate('TestStateDependentDelay', final_time = 0.75, maxh = 0.1)
250    t, x = res['time'], res['x']
251    x_expected = 1-t
252    x_expected2 = t-2+2*N.exp(0.5-t)
253    x_expected[t > 0.5] = x_expected2[t > 0.5]
254    assert_close(x, x_expected, 1e-3)
255
256@testattr(stddist_full = True)
257def test_delay_starting_at_zero():
258    res = compile_and_simulate('TestDelayStartingAtZero', final_time = 5, maxh = 0.1)
259    t, x = res['time'], res['x']
260    x_expected = N.exp(-t)
261    assert_close(x, x_expected, 1e-2)
262
263@testattr(stddist_full = True)
264def test_delay_starting_at_zero_no_event():
265    # Check that we get the same answer as above with noEvent
266    res = compile_and_simulate('TestDelayStartingAtZeroNoEvent', final_time = 5, maxh = 0.1)
267    t, x = res['time'], res['x']
268    x_expected = N.exp(-t)
269    assert_close(x, x_expected, 1e-2)
270
271@testattr(stddist_full = True)
272def test_variable_delay_events():
273    res = compile_and_simulate('TestVariableDelayEvents', final_time = 4, maxh = 0.1)
274    t, x = res['time'], res['x']
275    (inds,) = N.nonzero(N.diff(t) == 0)
276    x_expected = t-(N.cos(5*t)+1) > 1
277    x_expected[inds]   = x_expected[inds-1]
278    x_expected[inds+1] = x_expected[inds+2]
279    print(x_expected+0)
280    print(N.asarray(x, dtype=int))
281    print(t)
282    assert_close(x, x_expected, 1e-12)
283    assert_close(t[inds], N.array([0.43439307, 0.92808171, 1.0, 1.47239509,
284                                   1.64366485, 2.24956004, 2.67833465]), 1e-7)
285
286@testattr(stddist_full = True)
287def test_delay_going_to_zero():
288    res = compile_and_simulate('TestDelayGoingToZero', final_time = 2, maxh = 0.1)
289    t, x = res['time'], res['x']
290    x_expected = 0.22981923*N.exp(1-t)   
291    assert_close(x[t >= 1], x_expected[t >= 1], 1e-3)
292
293def check_zeno_repeat(t, x, ownevents=True):
294    (sw,)=N.nonzero(N.diff(t) == 0)
295    (nosw,)=N.nonzero(N.diff(t) > 0)
296    assert all((x==0) | (x == 1))
297    assert all(x[nosw]==x[nosw+1])
298    if ownevents: assert all((x[sw]!=x[sw+1]) | (N.abs(t[sw]-1)<1e-7))
299    tau = 1-t[sw]
300    if ownevents: assert_close(tau[1:], N.sqrt(0.5)*tau[0:-1], 1e-7)   
301
302@testattr(stddist_full = True)
303def test_zeno_repeat():
304    res = compile_and_simulate('TestZenoRepeat', final_time = 2, maxh = 0.1)
305    t, x = res['time'], res['x']
306    check_zeno_repeat(t, x)
307
308def check_zeno_repeat_noevent(t, x, ownevents=True):
309    points = t < 0.78
310    t, x = t[points], x[points]
311    xe1 = N.mod(N.log2(1-(t-1e-8)),1) > 0.5;
312    xe2 = N.mod(N.log2(1-(t+1e-8)),1) > 0.5;
313    (inds,) = N.nonzero(N.diff(t) == 0)
314
315    if ownevents: assert_close(t[inds], [1-N.sqrt(0.5)], 1e-7)
316    assert all(N.minimum(xe1, xe2) <= x)
317    assert all(x <= N.maximum(xe1, xe2))
318
319@testattr(stddist_full = True)
320def test_zeno_repeat_noevent():
321    res = compile_and_simulate('TestZenoRepeatNoEvent', final_time = 0.75, maxh = 1/20.5)
322    t, x = res['time'], res['x']
323    check_zeno_repeat_noevent(t, x)
324
325def sol_repeating_events2(t, d=1):
326    x_expected = N.mod(t, d)
327    (inds,) = N.nonzero(N.diff(t) == 0)
328    x_expected[inds]   = N.mod(t[inds]-1e-7, d)
329    x_expected[inds+1] = N.mod(t[inds]+1e-7, d)
330    return x_expected
331
332#@testattr(stddist_full = True)
333@nottest
334def test_multiple_delays():
335    res = compile_and_simulate('TestMultipleDelays', final_time = 10)
336    t   = res['time']
337    xr  = res['rep.x']
338    xrn = res['rep_ne.x']
339    xz  = res['zeno.x']
340    xzn = res['zeno_ne.x']
341
342    phi = (1 + N.sqrt(5))/2
343
344    assert sum(N.abs(xr-sol_repeating_events2(t)) > 1e-6) <= 8
345    assert sum(N.abs(xrn-sol_repeating_events2(t, phi)) > 0.2) <= 30
346    check_zeno_repeat(t/5, xz, False)
347    check_zeno_repeat_noevent(t/(5*phi), xzn, False)
348
349
350class TestSpatialDistribution:
351    def get_class_postfix(self):
352        return ""
353
354    def compile_and_load(self, class_name):
355        return compile_and_load(class_name+self.get_class_postfix(), filename=path_to_spatialdist_mo)
356
357    def compile_and_simulate(self, class_name, *args, **kwargs):
358        return compile_and_simulate(class_name+self.get_class_postfix(), *args, filename=path_to_spatialdist_mo, **kwargs)
359
360    @testattr(stddist_full = True)
361    def test_forward_flow(self):
362        res = self.compile_and_simulate('TestForwardFlow', final_time = 2, maxh = 0.01)
363        t, x = res['time'], res['x']
364        x_expected = switch_signal(t, N.sqrt(N.maximum(0, t**2-1)), 1+t*0, t-1);
365        assert_close(x, x_expected, 1e-2)
366
367    @testattr(stddist_full = True)
368    def test_back_flow(self):
369        res = self.compile_and_simulate('TestBackFlow', final_time = 3, maxh = 0.01)
370        t, x = res['time'], res['x']
371        x_expected = (t < 1)*t + (t >= 1)*((t<2)*(1-N.sqrt(N.maximum(0,t-1))) + (t>=2)*(2-t))
372        assert_close(x, x_expected, 1e-2)
373
374    @testattr(stddist_full = True)
375    def test_initial_contents(self):
376        res = self.compile_and_simulate('TestInitialContents', final_time = 2, maxh = 0.01)
377        t, x = res['time'], res['x']
378        x_expected = switch_signal(t, switch_signal(t, N.minimum(4*t, (0.75-t)/0.5), (1-t)/0.25, 0.75-t), -1, 1-t)
379        assert_close(x, x_expected, 1e-8)
380
381    @testattr(stddist_full = True)
382    def test_reversing_flow(self):
383        res = self.compile_and_simulate('TestReversingFlow', final_time = 10, maxh = 0.1)
384        t, x = res['time'], res['x']
385        inds = N.flatnonzero(N.diff(t) == 0)
386        x_expected = (N.cumsum(N.hstack([False, N.diff(t)==0])) & 3) > 0
387        assert_close(x, x_expected, 1e-15)
388        assert_close(t[inds], N.pi*N.array([1.0/2, 5.0/6, 1+1.0/2, 1+5.0/6, 2+1.0/2, 2+5.0/6]), 1e-7)
389
390    @testattr(stddist_full = True)
391    def test_sinusoid(self):
392        res = self.compile_and_simulate('TestSinusoid', final_time = 5, maxh = 0.1)
393        t, x = res['time'], res['x']
394        x_expected = 1.03*N.cos((t**2+0.35)*N.pi/2)
395        inds = (t >= 0.5)
396        assert_close(x[inds], x_expected[inds], 0.1)
397
398    @testattr(stddist_full = True)
399    def test_feed_loop(self, model_name = 'TestFeedLoop'):
400        res = self.compile_and_simulate(model_name, final_time = 10, maxh = 0.1)
401        t, x = res['time'], res['x']
402        x_expected = switch_signal(t, N.mod(4.5*N.sin(t)+0.5,2)-0.5, -(N.mod(4.5*N.sin(t)-0.5,2)-0.5), N.sin(N.pi*4.5*N.sin(t)))
403        assert_close(x, x_expected, 1e-8)
404
405    @testattr(stddist_full = True)
406    def test_feed_loop_no_pvel_events(self):
407        if self.get_class_postfix() != "":
408            def get_class_postfix():
409                return "(redeclare block SD=TestFeedLoopNoPVelEvents.SpatialDistReverse)"
410            self.get_class_postfix = get_class_postfix
411        self.test_feed_loop(model_name = 'TestFeedLoopNoPVelEvents')
412
413class TestSpatialDistributionRev(TestSpatialDistribution):
414    def get_class_postfix(self):
415        return "(redeclare block SD=SpatialDistReverse)"
Note: See TracBrowser for help on using the repository browser.