I am using scipy.solve_ivp to solve an ODE in a loop, such that in every iteration I change some term in the ODE (which is like the force applied).
In terms of my problem - I just cant seem to understand why when I track events of zero crossing the solution sometimes stalls and essentially zeros out (as the figure below shows). The code is a bit long, but only a small portion of it is actually relevant to my problem
(1) initialize some constants
MASS = 2.6e-4
WING_LENGTH = 0.07 # meters
GYRATION_RADIUS = 0.6 * WING_LENGTH
MoI = MASS * GYRATION_RADIUS ** 2 # Moment of Inertia
AIR_DENSITY = 1.2
WING_AREA = 0.5 * WING_LENGTH * (0.5 * WING_LENGTH) * np.pi # 1/2 ellipse with minor radios ~ 1/2 major = length/2
# drag coefficients from whitney & wood (JFM 2010):
C_D_MAX = 3.4
C_D_0 = 0.4
C_L_MAX = 1.8
ZERO_CROSSING = 1
RADIAN45 = np.pi / 4
RADIAN135 = 3 * RADIAN45
(2) Define a wrapper which is a RobotSimulation class - it's an easy interface that allows to easily change the values of the ODE. The important function is solve_dynamics (the others are not part of the problem I believe)
def phi_dot_zero_crossing_event(t, y):
"""
this event is given to solve_ivp to track if phi_dot == 0
:param t: unused time variable
:param y: a vector of [phi,phi_dot]
:return:
"""
return y[1]
class RobotSimulation:
def __init__(self, motor_torque=lambda x: 0, alpha=RADIAN45, phi0=0.0, phi_dot0=0.01, start_t=0, end_t=0.05, delta_t=0.001):
self.solution = 0
self.motor_torque = motor_torque
self.alpha = alpha # angle of attack
self.phi0 = phi0
self.phi_dot0 = phi_dot0
self.start_t = start_t
self.end_t = end_t
self.delta_t = delta_t
def set_time(self, new_start, new_end) -> None:
"""
we update the start time to be the end of the prev time and the end time
as the last end time + the window size we wish
"""
self.start_t = new_start
self.end_t = new_end
def set_init_cond(self, new_phi0, new_phi_dot0) -> None:
"""
we update the initial conditions to be the last value of previous iterations
"""
self.phi0 = new_phi0
self.phi_dot0 = new_phi_dot0
def set_motor_torque(self, new_motor_torque) -> None:
"""
sets the moment to the new value
"""
self.motor_torque = new_motor_torque
def flip_alpha(self) -> None:
"""
(Irrelevant for my question) changes the angle of attack's sign
:return:
"""
self.alpha = RADIAN135 if self.alpha == RADIAN45 else RADIAN45
def c_drag(self) -> float:
"""
calculates the drag coefficient based on the angle of attack
"""
return (C_D_MAX + C_D_0) / 2 - (C_D_MAX - C_D_0) / 2 * np.cos(2 * self.alpha)
def drag_torque(self, phi_dot):
"""
the drag moment
"""
return 0.5 * AIR_DENSITY * WING_AREA * self.c_drag() * GYRATION_RADIUS * phi_dot * np.abs(phi_dot)
def phi_derivatives(self, t, y):
"""
A function that defines the ODE that is to be solved: I * phi_ddot = tau_z - tau_drag.
We think of y as a vector y = [phi,phi_dot]. the ode solves dy/dt = f(y,t)
:return:
"""
phi, phi_dot = y[0], y[1]
dy_dt = [phi_dot, (self.motor_torque(t) - self.drag_torque(phi_dot)) / MoI]
return dy_dt
def solve_dynamics(self, phiarr, phidotarr, phiddotarr, angarr, torquearr):
"""
solves the ODE
:return:
"""
phi_0, phi_dot_0 = self.phi0, self.phi_dot0
start_t, end_t, delta_t = self.start_t, self.end_t, self.delta_t
phi_dot_zero_crossing_event.terminal = True
ang = []
times_between_zero_cross = []
sol_between_zero_cross = []
while start_t < end_t:
torquearr.append(self.motor_torque(start_t))
sol = solve_ivp(self.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0,phi_dot_0],events=phi_dot_zero_crossing_event)
self.solution = sol.y
ang.append(self.alpha * np.ones(len(sol.t))) # set alpha for every t based on solution's size
times_between_zero_cross.append(sol.t)
sol_between_zero_cross.append(sol.y)
if sol.status == ZERO_CROSSING:
start_t = sol.t[-1] + delta_t
phi_0, phi_dot_0 = sol.y[0][-1], sol.y[1][-1] # last step is now initial value
self.flip_alpha()
else: # no zero crossing
break
time = np.concatenate(times_between_zero_cross)
phi, phi_dot = np.concatenate(sol_between_zero_cross, axis=1)
ang = np.concatenate(ang)
_, phi_ddot = self.phi_derivatives(time, [phi, phi_dot])
phiarr.append(phi)
phidotarr.append(phi_dot)
phiddotarr.append(phi_ddot)
angarr.append(ang)
(3) Finally the main function calls solve_dynamics sequentially for different values of actions and updates the time window and initial conditions:
phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr = [], [], [], [], []
phi0, phidot0 = 0, 0.01
start_t, end_t, delta_t = 0, 0.05, 0.001
sim = RobotSimulation()
for action in [1, -1, 1, -1, 1, -1]:
sim.set_motor_torque(lambda x: action)
sim.solve_dynamics(phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr)
phi0, phidot0 = sim.solution[0][-1], sim.solution[1][-1]
start_t = end_t
end_t += 0.05
sim.set_init_cond(phi0, phidot0)
sim.set_time(start_t, end_t)
phi_arr = np.concatenate(phi_arr)
phi_dot_arr = np.concatenate(phi_dot_arr)
phi_ddot_arr = np.concatenate(phi_ddot_arr)
angle_arr = np.concatenate(angarr)
Intuitively, as the actions I apply are just repeating +1 and -1, I cannot understand why the blue points in the figure (see below) act so different
