8강 Bouncing ball, a hybrid system

Hybrid system


이전 시간, 라그랑지안을 통해 운동방정식을 세우는 방법에 대해서 학습하였는데요. 모든 로봇 시스템의 한 세트의 운동방정식으로 정의되면 좋겠지만, 벽과 충돌하거나, 강한 외력이 작용하면 시스템의 운동 상태가 변화하게 됩니다.

이를 로보틱스에서는 Hybrid system이라고 부르며, 바닥과 상호작용하는 보행 로봇, 형태가 변화하는 로봇 등에 사용됩니다.

Hybrid system의 학습을 위해 이번 시간 구현해볼 것은 바닥을 튀는 공입니다.

바닥과 충돌한 공은 탄성 계수만큼의 속도 감소 후 다시 위로 튀어오르고, 모든 에너지가 사라질때까지 이 과정을 반복합니다.

물체에 작용하는 힘은 중력과 더불어 속도의 제곱에 비례하는 항력(drag force)이 있으며, 이에 대한 운동방정식은 일전 계산한 바 있습니다. 따라서 이번에는 시스템의 상태가 변화하는 부분에 집중하고 구현해보겠습니다.

⇒ 상수 e는 물체와 바닥면과의 충돌 시 사용되는 탄성 계수입니다.

코드 구현

  • 구현할 코드를 우선 실행해볼까요?
python3 ball_bounce.py

⇒ 바닥과 공이 접촉하는 순간 시스템의 상태를 변화시켜야 합니다. 이러한 작업을 python에서는 어떻게 구현할 수 있는지 먼저 살펴보겠습니다.

solve ivp

python에서 미분 방정식의 해를 수치적으로 구하는 방법 중 하나인 solve_ivp에 대해 배워보겠습니다.

solve_ivp는 “initial value problem for a system of ODEs”의 약자로, 미분 방정식과 초기 조건이 주어지면 해당 조건을 바탕으로 수치적인 계산값을 도출해줍니다.

예시를 통해 python의 solve_ivp를 살펴봅시다.

  • 아래와 같은 미분방정식이 있을 때, solve_ivp를 통해 그래프를 그려보겠습니다.
# differential equation
def dydt(t, y):
    return -y + 2*np.sin(t)

전체 코드는 다음과 같습니다. 미분방정식을 구현한 함수, 계산이 적용되는 시간 범위와 초기값을 전달해주면 시간과 상태를 담은 solution을 반환합니다.


  • solve_ivp 매개변수

Parameter Description
f function
t_span 시간 범위
y0 초기값 (단일 변수일 시, list로 전달해야 합니다.
t_eval evaluation 할 범위로, 이 값은 적어도 t_span에 포함되는 범위여야 합니다.

  • 결과값

Parameter Description
sol.t 시간 범위
sol.y 계산 결과는 항상 list of list 입니다. 현재 변수가 하나뿐이므로 [0]으로 내부값을 가져왔습니다.


  • 이번에는 다변수함수의 미분방정식을 계산해봅시다.
def f(t, r):
    x, y = r
    
    fx = np.cos(y)
    fy = np.sin(x)
    
    return fx, fy


  • 다변수함수의 경우 initial state에 list가 아닌 tuple이 적용되며, 이는 방정식 함수 내부에서 따로 파싱됩니다.


solve_ivp - Event

앞으로 hybrid system의 운동방정식에 자주 사용될 기능으로, 특정 event가 발생하는 해를 계산하는 기능입니다.
  • 아래와 같이 y = 1을 만족하는 해를 event로 두고 이를 solve_ivp에게 전달해보겠습니다.
def event(t, y):
    return y[0] - 1

⇒ 미분방정식 함수와 달리 변수가 list 타입으로 들어오기 때문에 indexing을 해주어야 합니다.

  • y=1을 만족하는 해를 구하고 이를 시각화해봅시다. 전체 코드는 다음과 같습니다.

⇒ 지금은 특정 해를 구하는 상황이지만, 이후 예시들에서는 운동방정식의 조건이 바뀌는 event에 사용될 것입니다.

  • 결과 plot은 다음과 같습니다. span 범위 내에서 y = 1을 만족하는 지점들이 모두 계산되었습니다.


  • event condition - event에 terminal, direction을 부여하면 특정 조건을 만족하는 해를 분류할 수 있습니다.
def event(t, y):
    return y[0] - 1

event.terminal = False
event.direction = -1

Description
terminal 이벤트가 발생할 시 계산을 종료할지 여부입니다. 기본값은 False입니다.
direction 이 값이 +이면 이벤트 조건값 기준 음수에서 양수로 전환될 때만 terminal이 트리거되고, 반대로 방향이 -이면 이벤트 조건값 기준 양수에서 음수로 전환될 때 이벤트가 트리거됩니다. 0인 경우 방향에 상관 없이 terminal을 트리거합니다. 기본값은 0입니다.

⇒ 지금은 direction = -1이기 때문에 값이 감소하는 지점에서의 해만 추려낸 것이지요.


  • 다변수 함수의 event도 구현해봅시다.

terminal 조건은 따로 명시하지 않았기 때문에 event를 만족하는 모든 해가 표시되었고, 범위 0 ~ 10 사이 y=1을 만족하는 해가 모두 표시되었습니다.


Bouncing Ball 코드 구현

코드의 전체 흐름은 아래와 같은 맥락으로 흘러갑니다.

  1. simulation
    1. one_bounce
      1. solve_ivp
        1. projectile - 공의 운동방정식을 포함하고 있습니다.
      2. contact - 공과 바닥이 접촉하는 event입니다.
  2. animate
  3. plot


  • projectile - 운동방정식이 구현되어 있는 함수로 $x, dx, y, dy$를 input으로 받아 ⇒ $dx, ddx, dy, ddy$를 output으로 리턴합니다.


  • contact - solve_ivp의 event 함수로 바닥에 공이 닿는 시점인 y = 0이 return 값이 됩니다.


  • solve_ivp - 실제 미분방정식의 수치적 계산이 이루어지며, 매번 bounce가 반복될 때마다 공의 상태 벡터를 계산해줄 것입니다.

⇒ 여기에 사용된 추가 매개변수들에 대해서 좀 더 알고 싶다면 solve_ivp의 공식 문서를 살펴보시기 바랍니다. 최적화에 대해 보다 깊은 이해를 가질 수 있습니다.


  • one_bounce - 한번의 bounce를 함수화해둔 구현입니다.

이 함수에서 contact event에 대한 조건이 부여됩니다. contact.terminal을 통해 event 발생 시 계산이 멈출지 여부를 정하고 , contact.direction이 -1이기 때문에 event가 발생하는 조건은 y값이 감소할 때가 됩니다.

더불어, contact이 발생하면 y축 속도의 방향과 크기가 변경되기 때문에 초기 조건 z0를 다시 설정하고 있습니다.


  • 지금까지 살펴본 코드가 one_bounce였습니다. 이제 특정 시간동안 bounce가 반복되는 상황을 시뮬레이션 해봅시다. one_bounce가 반복될 때마다 해당 값들을 차곡차곡 배열에 저장하며,


  • 이렇게 누적된 시뮬레이션 데이터들은 이후 animation과 Plotting에 사용됩니다.
t, z = simulation(t0, t_end, z0, params)


이번 시간 살펴본 hybrid 시스템은 이후 legged robotics 강의에서 매번 사용되므로 두 번 이상 복습을 추천드립니다!

Complete and Continue