(Q4)The following parameters are given: M = 0.3 kg, (5.3a) m = 0.1 kg, (5.3b) ` = 35 cm, (5.3c) g = 9.81 m/s 2 . (5.3d) Suppose that the inverted pendulum is initially at rest at the upright position. Simulate the trajectories of θ(t) and x(t) upon the action of the force F(t) = sin(100t 2 ). (5.4) Simulate over t ∈ [0, 0.2]s.

PendulumQuestion4.py

import sympy as sym
import numpy as np
import control as ctrl
import matplotlib.pyplot as p
m, ell, x3, x4, M, g, F, m = sym.symbols('m, ell, x3, x4, M, g, F, m')

phi = 4*m*ell*x4**2*sym.sin(x3) + 4*F - 3*m*g*sym.sin(x3)*sym.cos(x3)
phi /= 4*(M+m) - 3*m*sym.cos(x3)**2
dphi_x3 = phi.diff(x3)
dphi_x4 = phi.diff(x4)
dphi_F = phi.diff(F)

psi = -3*(m*ell*x4**2*sym.sin(x3)*sym.cos(x3) + F*sym.cos(x3) - (M+m)*g*sym.sin(x3))
psi /= (4*(M+m) - 3*m*sym.cos(x3)**2)*ell
dpsi_x3 = psi.diff(x3)
dpsi_x4 = psi.diff(x4)
dpsi_F = psi.diff(F)
# Equilibrium point
Feq = 0
x3eq = 0
x4eq = 0
#------------------------------
#  adding equilibrium to x1 and x2
#------------------------------
x1eq = 0
x2eq = 0

dphi_F_eq = dphi_F.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dphi_x3_eq = dphi_x3.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dphi_x4_eq = dphi_x4.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])

dpsi_F_eq = dpsi_F.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dpsi_x3_eq = dpsi_x3.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dpsi_x4_eq = dpsi_x4.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])

a = dphi_F_eq
b = -dphi_x3_eq
c = 3/(ell*(4*M + m))
d = 3*(M+m)*g/(ell*(4*M + m))
# GIVEN VALUES!


def evaluate_at_given_parameters(z):

    M_value = 0.3
    m_value = 0.1
    g_value = 9.81
    ell_value = 0.35

    return float(z.subs([(M, M_value), (m, m_value), (ell, ell_value), (g, g_value)]))

a_value = evaluate_at_given_parameters(a)
b_value = evaluate_at_given_parameters(b)
c_value = evaluate_at_given_parameters(c)
d_value = evaluate_at_given_parameters(d)

# -----------------------------------
# Control library of Python
# -----------------------------------
a, b, c, d = sym.symbols('a:d', real=True, positive=True)
s, t = sym.symbols('s, t')

transfer_function_F_to_x3 = -c/(s**2 - d)
transfer_function_F_to_x1 = (a*(s**2 - d) + b*c)/(s**2)*(s**2-d)

n_points = 500
t_final = 0.2
t_span = np.linspace(0, t_final, n_points)
input_signal = np.sin(100*(t_span**2))
transfer_function_F2x3 = ctrl.TransferFunction([-c_value], [1, 0, -d_value])
transfer_function_F2x1 = ctrl.TransferFunction([a_value, 0, (-a_value*d_value) + (b_value*c_value)], [1, 0, -d_value, 0, 0])

tf_out1, t_out1, Ft_out1= ctrl.forced_response(transfer_function_F2x3, t_span, input_signal)
tf_out2, t_out2, Ft_out2= ctrl.forced_response(transfer_function_F2x1, t_span, input_signal)

plt.plot(tf_out1, t_out1)
plt.xlabel('Time (s)')
plt.ylabel('Angle (rad)')
plt.grid()
plt.savefig("Question4Angle(rads)AgainstTime.eps", format = "eps")
plt.show()
plt.plot(tf_out2, t_out2)
plt.xlabel('Time (s)')
plt.ylabel('Horizontal Position (m)')
plt.grid()
plt.savefig("Question4HorizontalPosistionAgainstTime.eps", format = "eps")
plt.show()