-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathdampened_oscillator_declarative.py
130 lines (91 loc) · 3.52 KB
/
dampened_oscillator_declarative.py
1
2
3
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
from numerous.declarative.specification import ScopeSpec, ItemsSpec, Module, EquationSpec
from numerous.declarative.variables import Parameter, Constant, State
from numerous.declarative.bus import Connector, create_connections, get_value_for, set_value_from
import numpy as np
class DampenedOscillator(Module):
"""
Equation and item modelling a spring and dampener
"""
tag = "dampened_oscillator"
class Mechanics(ScopeSpec):
# Define variables
k = Constant(1) # spring constant
c = Constant(1) # coefficient of friction
a = Parameter(0) # acceleration
x, x_dot = State(0) # distance
v, v_dot = State(0) # velocity
mechanics = Mechanics()
coupling = Connector(
x = set_value_from(mechanics.x),
F = get_value_for(mechanics.v_dot)
)
@EquationSpec(mechanics)
def eval(self, scope: Mechanics):
scope.a = -scope.k * scope.x - scope.c * scope.v
scope.v_dot = scope.a
scope.x_dot = scope.v
class SpringCoupling(Module):
tag = "springcoup"
class Mechanics(ScopeSpec):
k = Parameter(1)
c = Parameter(1)
F1 = Parameter(0)
F2 = Parameter(0)
x1 = Parameter(0)
x2 = Parameter(0)
mechanics = Mechanics()
class Items(ItemsSpec):
side1: DampenedOscillator
side2: DampenedOscillator
side1 = Connector(
F = set_value_from(mechanics.F1),
x = get_value_for(mechanics.x1)
)
side2 = Connector(
F=set_value_from(mechanics.F2),
x=get_value_for(mechanics.x2)
)
items = Items()
def __init__(self, tag="springcoup", k=1.0):
super().__init__(tag)
self.mechanics.set_values(k=k)
@EquationSpec(mechanics)
def eval(self, scope: Mechanics):
scope.c = scope.k
dx = scope.x1 - scope.x2
F = np.abs(dx) * scope.c
scope.F1 = -F if scope.x1 > scope.x2 else F # [kg/s]
scope.F2 = -scope.F1
class OscillatorSystem(Module):
class Items(ItemsSpec):
# The oscillator system will have two oscillators interacting via the coupling
oscillator1: DampenedOscillator
oscillator2: DampenedOscillator
coupling: SpringCoupling
items = Items()
with create_connections() as connections:
# connect oscillators via the coupling
items.oscillator1.coupling >> items.coupling.side1
items.oscillator2.coupling >> items.coupling.side2
def __init__(self, k=0.01, c=0.001, x0=(1, 2), tag=None):
super(OscillatorSystem, self).__init__(tag=tag)
# Initialized the modules in the oscillator system
self.items.oscillator1 = DampenedOscillator(tag='osc1')
self.items.oscillator1.mechanics.set_values(k=k, c=c, x=x0[0])
self.items.oscillator2 = DampenedOscillator(tag='osc2')
self.items.oscillator2.mechanics.set_values(k=k, c=c, x=x0[1])
self.items.coupling = SpringCoupling(k=k)
if __name__ == "__main__":
from numerous.engine import simulation
from numerous.engine import model
from matplotlib import pyplot as plt
subsystem = OscillatorSystem(tag='system', k=0.01, c=0.001, x0=[1.0, 2.0])
# Define simulation
s = simulation.Simulation(model.Model(subsystem, use_llvm=False), t_start=0, t_stop=500.0, num=1000, num_inner=100,
max_step=1)
# Solve and plot
s.solve()
s.model.historian_df[['system.osc1.mechanics.x', 'system.osc2.mechanics.x']].plot()
# print()
plt.show()
plt.interactive(False)