BM_15.py 5.38 KB
Newer Older
navidmokh's avatar
navidmokh committed
1
2
from scipy.integrate import odeint
import numpy as np
3
# import matplotlib.pyplot as plt
navidmokh's avatar
navidmokh committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

# Source: https://ths.rwth-aachen.de/research/projects/hypro/brusselator/

def BM15_dynamic(y, t):
    x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15 = y
    
    u1 = 0.9
    x1_dot = 0.0061322*u1 - 0.0075102*x1 + 5.2756*x2 - 0.00096238*x3 - 0.6634*x4 - 0.006689*x5 + 0.16081*x6 + 0.070705*x7 + 0.01843*x8 + 0.033143*x9 + 0.025063*x10 + 0.1069*x11 - 0.0095973*x12 + 0.054215*x13 + 0.0078958*x14 - 0.05736*x15
    x2_dot = 0.064531*u1 - 5.2756*x1 - 0.85737*x2 + 0.090377*x3 + 0.92166*x4 + 0.13158*x5 - 0.97335*x6 - 0.44973*x7 - 0.33292*x8 - 0.25467*x9 - 0.38152*x10 - 0.92917*x11 + 0.12951*x12 - 0.47949*x13 - 0.10417*x14 + 0.51986*x15
    x3_dot = 0.00069603*u1 - 0.00096238*x1 - 0.090377*x2 - 0.00012541*x3 - 13.542*x4 - 0.00092367*x5 + 0.026749*x6 + 0.011595*x7 + 0.0025622*x8 + 0.0051874*x9 + 0.0035389*x10 + 0.016317*x11 - 0.0013727*x12 + 0.0082485*x13 + 0.0011323*x14 - 0.0086882*x15
    x4_dot = 0.6634*x1 - 0.062226*u1 + 0.92166*x2 + 13.542*x3 - 1.004*x4 - 0.17895*x5 + 1.117*x6 + 0.51795*x7 + 0.44311*x8 + 0.29783*x9 + 0.4893*x10 + 1.0996*x11 - 0.1625*x12 + 0.56846*x13 + 0.13022*x14 - 0.61793*x15
    x5_dot = 0.0035048*u1 - 0.006689*x1 - 0.13158*x2 - 0.00092367*x3 + 0.17895*x4 - 0.008656*x5 + 23.761*x6 + 1.1717*x7 + 0.024823*x8 + 0.13715*x9 + 0.037225*x10 + 0.3119*x11 - 0.015556*x12 + 0.15229*x13 + 0.013037*x14 - 0.15325*x15
    x6_dot = 0.047228*u1 - 0.16081*x1 - 0.97335*x2 - 0.026749*x3 + 1.117*x4 - 23.761*x5 - 1.5873*x6 - 0.75053*x7 - 4.9343*x8 - 0.4724*x9 - 1.9409*x10 - 1.8836*x11 + 0.4708*x12 - 0.98584*x13 - 0.36077*x14 + 1.0912*x15
    x7_dot = 0.021423*u1 - 0.070705*x1 - 0.44973*x2 - 0.011595*x3 + 0.51795*x4 - 1.1717*x5 - 0.75053*x6 - 0.35552*x7 - 6.1681*x8 - 0.22575*x9 - 1.0964*x10 - 0.90761*x11 + 0.24584*x12 - 0.47568*x13 - 0.18672*x14 + 0.52758*x15
    x8_dot = 0.01843*x1 - 0.0093666*u1 + 0.33292*x2 + 0.0025622*x3 - 0.44311*x4 + 0.024823*x5 + 4.9343*x6 + 6.1681*x7 - 0.071566*x8 - 0.55207*x9 - 0.10886*x10 - 1.0701*x11 + 0.046135*x12 - 0.51571*x13 - 0.038788*x14 + 0.5106*x15
    x9_dot = 0.01125*u1 - 0.033143*x1 - 0.25467*x2 - 0.0051874*x3 + 0.29783*x4 - 0.13715*x5 - 0.4724*x6 - 0.22575*x7 + 0.55207*x8 - 0.14992*x9 - 14.126*x10 - 0.62987*x11 + 0.31572*x12 - 0.33265*x13 - 0.21843*x14 + 0.37322*x15
    x10_dot = 0.025063*x1 - 0.011919*u1 + 0.38152*x2 + 0.0035389*x3 - 0.4893*x4 + 0.037225*x5 + 1.9409*x6 + 1.0964*x7 - 0.10886*x8 + 14.126*x9 - 0.17215*x10 - 3.3095*x11 + 0.075998*x12 - 1.4714*x13 - 0.064505*x14 + 1.3304*x15
    x11_dot = 0.038892*u1 - 0.1069*x1 - 0.92917*x2 - 0.016317*x3 + 1.0996*x4 - 0.3119*x5 - 1.8836*x6 - 0.90761*x7 + 1.0701*x8 - 0.62987*x9 + 3.3095*x10 - 2.775*x11 + 33.887*x12 - 1.4783*x13 - 4.2421*x14 + 1.681*x15
    x12_dot = 0.0043369*u1 - 0.0095973*x1 - 0.12951*x2 - 0.0013727*x3 + 0.1625*x4 - 0.015556*x5 - 0.4708*x6 - 0.24584*x7 + 0.046135*x8 - 0.31572*x9 + 0.075998*x10 - 33.887*x11 - 0.035149*x12 + 5.3142*x13 + 0.030178*x14 - 2.004*x15
    x13_dot = 0.019912*u1 - 0.054215*x1 - 0.47949*x2 - 0.0082485*x3 + 0.56846*x4 - 0.15229*x5 - 0.98584*x6 - 0.47568*x7 + 0.51571*x8 - 0.33265*x9 + 1.4714*x10 - 1.4783*x11 - 5.3142*x12 - 0.78885*x13 - 6.3028*x14 + 0.8993*x15
    x14_dot = 0.0078958*x1 - 0.0035327*u1 + 0.10417*x2 + 0.0011323*x3 - 0.13022*x4 + 0.013037*x5 + 0.36077*x6 + 0.18672*x7 - 0.038788*x8 + 0.21843*x9 - 0.064505*x10 + 4.2421*x11 + 0.030178*x12 + 6.3028*x13 - 0.025985*x14 + 4.0141*x15
    x15_dot = 0.05736*x1 - 0.021347*u1 + 0.51986*x2 + 0.0086882*x3 - 0.61793*x4 + 0.15325*x5 + 1.0912*x6 + 0.52758*x7 - 0.5106*x8 + 0.37322*x9 - 1.3304*x10 + 1.681*x11 + 2.004*x12 + 0.8993*x13 - 4.0141*x14 - 1.0293*x15

    dydt = [x1_dot,x2_dot,x3_dot,x4_dot,x5_dot,x6_dot,x7_dot,x8_dot,x9_dot,x10_dot,x11_dot,x12_dot,x13_dot,x14_dot,x15_dot]
    return dydt

def TC_Simulate(Mode,initialCondition,time_bound):
    time_step = 0.005;
    time_bound = float(time_bound)

    number_points = int(np.ceil(time_bound/time_step))
    t = [i*time_step for i in range(0,number_points)]
    if t[-1] != time_step:
        t.append(time_bound)
    newt = []
    for step in t:
        newt.append(float(format(step, '.3f')))
    t = newt

    sol = odeint(BM15_dynamic, initialCondition, t, hmax=time_step)

    # Construct the final output
    trace = []
    for j in range(len(t)):
        #print t[j], current_psi
        tmp = []
        tmp.append(t[j])
        tmp.append(float(sol[j,0]))
        tmp.append(float(sol[j,1]))
        tmp.append(float(sol[j,2]))
        tmp.append(float(sol[j,3]))
        tmp.append(float(sol[j,4]))
        tmp.append(float(sol[j,5]))
        tmp.append(float(sol[j,6]))
        tmp.append(float(sol[j,7]))
        tmp.append(float(sol[j,8]))
        tmp.append(float(sol[j,9]))
        tmp.append(float(sol[j,10]))
        tmp.append(float(sol[j,11]))
        tmp.append(float(sol[j,12]))
        tmp.append(float(sol[j,13]))
        tmp.append(float(sol[j,14]))
        trace.append(tmp)
    return trace

if __name__ == "__main__":

    sol = TC_Simulate("Default", [0.0060145, -0.0009624, 0.0177578, -0.0047553, 0.0177717, -0.0122932, -0.0026115, -0.0122392, 0.0014043, -0.0219818, -0.0396104, -0.0042632, -0.0070418, -0.009347, 0.0022889], 10.0)
    for s in sol:
li213's avatar
li213 committed
73
        print(s)
navidmokh's avatar
navidmokh committed
74
75
76
77
78
79
80
81
82
83
84
85

    # time = [row[0] for row in sol]

    # a = [row[1] for row in sol]

    # b = [row[2] for row in sol]

    # plt.plot(time, a, "-r")
    # plt.plot(time, b, "-g")
    # plt.show()
    # plt.plot(a, b, "-r")
    # plt.show()