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

Last change on this file since 13800 was 13800, checked in by randersson, 8 weeks ago

#5819 Merged trunk into branch

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