Skip to content

Commit d6635f8

Browse files
JianqiangDingJianqiangDing
committed
update examples
1 parent 5a1d1e3 commit d6635f8

6 files changed

Lines changed: 111 additions & 103 deletions

File tree

‎example/benchmark_brusselator_cmp.py‎

Lines changed: 33 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -5,45 +5,47 @@
55
from pybdr.algorithm import ASB2008CDC
66
from pybdr.dynamic_system import NonLinSys
77
from pybdr.geometry import Zonotope, Interval, Geometry
8-
from pybdr.geometry.operation import boundary
8+
from pybdr.geometry.operation import boundary, cvt2
99
from pybdr.model import *
1010
from pybdr.util.visualization import plot, plot_cmp
1111

12-
# init dynamic system
13-
system = NonLinSys(Model(brusselator, [2, 1]))
12+
if __name__ == '__main__':
1413

15-
# settings for the computation
16-
options = ASB2008CDC.Options()
17-
options.t_end = 6.0
18-
options.step = 0.02
19-
options.tensor_order = 2
20-
options.taylor_terms = 4
14+
# init dynamic system
15+
system = NonLinSys(Model(brusselator, [2, 1]))
2116

22-
options.u = Zonotope([0], np.diag([0]))
23-
options.u_trans = options.u.c
17+
# settings for the computation
18+
options = ASB2008CDC.Options()
19+
options.t_end = 6.0
20+
options.step = 0.02
21+
options.tensor_order = 2
22+
options.taylor_terms = 4
2423

25-
# settings for the using geometry
26-
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
27-
Zonotope.ORDER = 50
24+
options.u = Zonotope([0], np.diag([0]))
25+
options.u_trans = options.u.c
2826

29-
z = Zonotope([0.2, 0.2], np.diag([0.1, 0.1]))
27+
# settings for the using geometry
28+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
29+
Zonotope.ORDER = 50
3030

31-
no_boundary_analysis = True
31+
z = Interval.identity(2) * 0.1 + 0.2;
3232

33-
tp_whole, tp_bound = None, None
34-
xlim, ylim = [0, 2], [0, 2.4]
33+
no_boundary_analysis = False
3534

36-
if no_boundary_analysis:
37-
# reachable sets computation without boundary analysis
38-
options.r0 = [z]
39-
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
40-
else:
41-
# reachable sets computation with boundary analysis
42-
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
43-
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
35+
tp_whole, tp_bound = None, None
36+
xlim, ylim = [0, 2], [0, 2.4]
4437

45-
# visualize the results
46-
if no_boundary_analysis:
47-
plot(tp_whole, [0, 1], xlim=xlim, ylim=ylim)
48-
else:
49-
plot(tp_bound, [0, 1], xlim=xlim, ylim=ylim)
38+
if no_boundary_analysis:
39+
# reachable sets computation without boundary analysis
40+
options.r0 = [cvt2(z, Geometry.TYPE.ZONOTOPE)]
41+
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
42+
else:
43+
# reachable sets computation with boundary analysis
44+
options.r0 = boundary(z, 0.3, Geometry.TYPE.ZONOTOPE)
45+
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
46+
47+
# visualize the results
48+
if no_boundary_analysis:
49+
plot(tp_whole, [0, 1], xlim=xlim, ylim=ylim)
50+
else:
51+
plot(tp_bound, [0, 1], xlim=xlim, ylim=ylim)

‎example/benchmark_jet_engine_cmp.py‎

Lines changed: 24 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -5,36 +5,37 @@
55
from pybdr.algorithm import ASB2008CDC
66
from pybdr.dynamic_system import NonLinSys
77
from pybdr.geometry import Zonotope, Interval, Geometry
8-
from pybdr.geometry.operation import boundary
8+
from pybdr.geometry.operation import boundary, cvt2
99
from pybdr.model import *
1010
from pybdr.util.visualization import plot, plot_cmp
1111

12-
# init dynamic system
13-
system = NonLinSys(Model(jet_engine, [2, 1]))
12+
if __name__ == '__main__':
13+
# init dynamic system
14+
system = NonLinSys(Model(jet_engine, [2, 1]))
1415

15-
# settings for the computation
16-
options = ASB2008CDC.Options()
17-
options.t_end = 10
18-
options.step = 0.01
19-
options.tensor_order = 2
20-
options.taylor_terms = 4
16+
# settings for the computation
17+
options = ASB2008CDC.Options()
18+
options.t_end = 10
19+
options.step = 0.01
20+
options.tensor_order = 2
21+
options.taylor_terms = 4
2122

22-
options.u = Zonotope([0], np.diag([0]))
23-
options.u_trans = options.u.c
23+
options.u = Zonotope([0], np.diag([0]))
24+
options.u_trans = options.u.c
2425

25-
# settings for the using geometry
26-
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
27-
Zonotope.ORDER = 50
26+
# settings for the using geometry
27+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
28+
Zonotope.ORDER = 50
2829

29-
z = Zonotope([1, 1], np.diag([0.2, 0.2]))
30+
z = Interval.identity(2) * 0.2 + 1
3031

31-
# reachable sets computation without boundary analysis
32-
options.r0 = [z]
33-
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
32+
# reachable sets computation without boundary analysis
33+
options.r0 = [cvt2(z, Geometry.TYPE.ZONOTOPE)]
34+
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
3435

35-
# reachable sets computation with boundary analysis
36-
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
37-
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
36+
# reachable sets computation with boundary analysis
37+
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
38+
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
3839

39-
# visualize the results
40-
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])
40+
# visualize the results
41+
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])

‎example/benchmark_lotka_volterra_2d_cmp.py‎

Lines changed: 24 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -5,36 +5,37 @@
55
from pybdr.algorithm import ASB2008CDC
66
from pybdr.dynamic_system import NonLinSys
77
from pybdr.geometry import Zonotope, Interval, Geometry
8-
from pybdr.geometry.operation import boundary
8+
from pybdr.geometry.operation import boundary, cvt2
99
from pybdr.model import *
1010
from pybdr.util.visualization import plot, plot_cmp
1111

12-
# init dynamic system
13-
system = NonLinSys(Model(lotka_volterra_2d, [2, 1]))
12+
if __name__ == '__main__':
13+
# init dynamic system
14+
system = NonLinSys(Model(lotka_volterra_2d, [2, 1]))
1415

15-
# settings for the computation
16-
options = ASB2008CDC.Options()
17-
options.t_end = 2.2
18-
options.step = 0.005
19-
options.tensor_order = 3
20-
options.taylor_terms = 4
16+
# settings for the computation
17+
options = ASB2008CDC.Options()
18+
options.t_end = 2.2
19+
options.step = 0.005
20+
options.tensor_order = 3
21+
options.taylor_terms = 4
2122

22-
options.u = Zonotope.zero(1, 1)
23-
options.u_trans = np.zeros(1)
23+
options.u = Zonotope.zero(1, 1)
24+
options.u_trans = np.zeros(1)
2425

25-
# settings for the using geometry
26-
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
27-
Zonotope.ORDER = 50
26+
# settings for the using geometry
27+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
28+
Zonotope.ORDER = 50
2829

29-
z = Zonotope([3, 3], np.diag([0.5, 0.5]))
30+
z = Interval([2.5, 2.5], [3.5, 3.5])
3031

31-
# reachable sets computation without boundary analysis
32-
options.r0 = [z]
33-
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
32+
# reachable sets computation without boundary analysis
33+
options.r0 = [cvt2(z, Geometry.TYPE.ZONOTOPE)]
34+
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
3435

35-
# reachable sets computation with boundary analysis
36-
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
37-
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
36+
# reachable sets computation with boundary analysis
37+
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
38+
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
3839

39-
# visualize the results
40-
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])
40+
# visualize the results
41+
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])

‎example/benchmark_synchronous_machine_cmp.py‎

Lines changed: 25 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -5,36 +5,38 @@
55
from pybdr.algorithm import ASB2008CDC
66
from pybdr.dynamic_system import NonLinSys
77
from pybdr.geometry import Zonotope, Interval, Geometry
8-
from pybdr.geometry.operation import boundary
8+
from pybdr.geometry.operation import boundary, cvt2
99
from pybdr.model import *
1010
from pybdr.util.visualization import plot, plot_cmp
1111

12-
# init dynamic system
13-
system = NonLinSys(Model(synchronous_machine, [2, 1]))
12+
if __name__ == '__main__':
13+
# init dynamic system
14+
system = NonLinSys(Model(synchronous_machine, [2, 1]))
1415

15-
# settings for the computation
16-
options = ASB2008CDC.Options()
17-
options.t_end = 3
18-
options.step = 0.01
19-
options.tensor_order = 3
20-
options.taylor_terms = 4
16+
# settings for the computation
17+
options = ASB2008CDC.Options()
18+
options.t_end = 3
19+
options.step = 0.01
20+
options.tensor_order = 3
21+
options.taylor_terms = 4
2122

22-
options.u = Zonotope.zero(1, 1)
23-
options.u_trans = np.zeros(1)
23+
options.u = Zonotope.zero(1, 1)
24+
options.u_trans = np.zeros(1)
2425

25-
# settings for the using geometry
26-
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
27-
Zonotope.ORDER = 50
26+
# settings for the using geometry
27+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
28+
Zonotope.ORDER = 50
2829

29-
z = Zonotope([0, 3], np.diag([0.2, 0.2]))
30+
# z = Zonotope([0, 3], np.diag([0.2, 0.2]))
31+
z = Interval([-0.2, 2.8], [0.2, 3.2])
3032

31-
# reachable sets computation without boundary analysis
32-
options.r0 = [z]
33-
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
33+
# reachable sets computation without boundary analysis
34+
options.r0 = [cvt2(z, Geometry.TYPE.ZONOTOPE)]
35+
ti_whole, tp_whole, _, _ = ASB2008CDC.reach(system, options)
3436

35-
# reachable sets computation with boundary analysis
36-
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
37-
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
37+
# reachable sets computation with boundary analysis
38+
options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
39+
ti_bound, tp_bound, _, _ = ASB2008CDC.reach(system, options)
3840

39-
# visualize the results
40-
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])
41+
# visualize the results
42+
plot_cmp([tp_whole, tp_bound], [0, 1], cs=["#FF5722", "#303F9F"])

‎pybdr/algorithm/asb2008cdc.py‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
from dataclasses import dataclass
1212

1313
import numpy as np
14+
15+
np.seterr(divide='ignore', invalid='ignore')
1416
from scipy.special import factorial
1517

1618
from pybdr.dynamic_system import LinSys, NonLinSys
@@ -169,7 +171,7 @@ def reach_one_step(cls, sys: NonLinSys, r0, err, opt: Options):
169171
temp_r_ti, temp_r_tp, abst_err, dims = cls.linear_reach(
170172
sys, r0[i], err[i], opt
171173
)
172-
# check if initial set has to be split
174+
# check if the initial set has to be split
173175
if len(dims) <= 0:
174176
r_ti.append(temp_r_ti)
175177
r_tp.append(temp_r_tp)
@@ -190,7 +192,7 @@ def reach(cls, sys: NonLinSys, opt: Options):
190192
next_ti, next_tp = cls.reach_one_step(sys, tp_set[-1], err[-1], opt)
191193
opt.step_idx += 1
192194
ti_set.append(next_ti)
193-
ti_time.append(time_pts[opt.step_idx - 1 : opt.step_idx + 1])
195+
ti_time.append(time_pts[opt.step_idx - 1: opt.step_idx + 1])
194196
tp_set.append(next_tp)
195197
tp_time.append(time_pts[opt.step_idx])
196198

‎pybdr/util/visualization/plot.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -161,7 +161,7 @@ def __2d_plot_cmp(collections, dims, width, height, xlim, ylim, cs, filled):
161161
raise NotImplementedError
162162

163163
ax.autoscale_view()
164-
ax.axis("equal")
164+
# ax.axis("equal")
165165
ax.set_xlabel("x" + str(dims[0]))
166166
ax.set_ylabel("x" + str(dims[1]))
167167

0 commit comments

Comments
 (0)