nonlinear physical system, rigid pendulum in which the pivot point vibrates in a vertical direction
simulation.py
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 7 19:06:37 2021
@author: Rocky
"""
import numpy as np
import system_model
from scipy.integrate import solve_ivp
from ackrep_core import ResultContainer
from ackrep_core.system_model_management import save_plot_in_dir
import matplotlib.pyplot as plt
import os
def simulate():
model = system_model.Model()
rhs_xx_pp_symb = model.get_rhs_symbolic()
print("Computational Equations:\n")
for i, eq in enumerate(rhs_xx_pp_symb):
print(f"dot_x{i+1} =", eq)
rhs = model.get_rhs_func()
# Initial State values
xx0 = [200 / 360 * 2 * np.pi, 0]
t_end = 10
tt = times = np.linspace(0, t_end, 10000)
sim = solve_ivp(rhs, (0, t_end), xx0, t_eval=tt)
uu = model.uu_func(sim.t, xx0)[0] * 0.005 + 180
sim.uu = uu
save_plot(sim)
return sim
def save_plot(sim):
# create figure + 2x2 axes array
fig1, axs = plt.subplots(nrows=2, ncols=1, figsize=(12.8, 9.6))
# print in axes top left
axs[0].plot(sim.t, np.real(sim.y[0] * 360 / (2 * np.pi)), label=r"$\varphi$")
axs[0].plot(sim.t, list(sim.uu), label=r"periodic excitation $\cos(\omega*t)$")
axs[0].set_ylabel("Angle [rad]") # y-label
axs[0].grid()
axs[0].legend()
# print in axes top right
axs[1].plot(sim.t, np.real(sim.y[1]))
axs[1].set_ylabel("Angular velocity [rad/s]") # y-label
axs[1].set_xlabel("Time [s]") # x-Label
axs[1].grid()
plt.tight_layout()
## static
save_plot_in_dir()
def evaluate_simulation(simulation_data):
"""
:param simulation_data: simulation_data of system_model
:return:
"""
expected_final_state = [3.1432908256013783, -0.014961262100684384]
rc = ResultContainer(score=1.0)
simulated_final_state = simulation_data.y[:, -1]
rc.final_state_errors = [
simulated_final_state[i] - expected_final_state[i] for i in np.arange(0, len(simulated_final_state))
]
rc.success = np.allclose(expected_final_state, simulated_final_state, rtol=0, atol=1e-2)
return rc
system_model.py
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 9 13:33:34 2021
@author: Jonathan Rockstroh
"""
import sympy as sp
import symbtools as st
import importlib
import sys, os
from ipydex import IPS, activate_ips_on_exception # for debugging only
from ackrep_core.system_model_management import GenericModel, import_parameters
# Import parameter_file
params = import_parameters()
class Model(GenericModel):
def initialize(self):
"""
this function is called by the constructor of GenericModel
:return: None
"""
# Define number of inputs
self.u_dim = 1
# Set "sys_dim" to constant value, if system dimension is constant
self.sys_dim = 2
# check existance of params file -> if not: System is defined to hasn't
# parameters
self.has_params = True
self.params = params
# ----------- SET DEFAULT INPUT FUNCTION ---------- #
# --------------- Only for non-autonomous Systems
# --------------- MODEL DEPENDENT
def uu_default_func(self):
"""
:return:(function with 2 args - t, xx_nv) default input function
"""
a, omega = self.pp_symb[2], self.pp_symb[3]
u_sp = self.pp_dict[a] * sp.sin(self.pp_dict[omega] * self.t_symb - sp.pi / 2)
du_dtt_sp = u_sp.diff(self.t_symb, 2)
du_dtt_sp = du_dtt_sp.subs(self.pp_subs_list)
du_dtt_func = st.expr_to_func(self.t_symb, du_dtt_sp)
def uu_rhs(t, xx_nv):
du_dtt_nv = du_dtt_func(t)
return [du_dtt_nv]
return uu_rhs
# ----------- SYMBOLIC RHS FUNCTION ---------- #
# --------------- MODEL DEPENDENT
def get_rhs_symbolic(self):
"""
:return:(matrix) symbolic rhs-functions
"""
if self.dxx_dt_symb is not None:
return self.dxx_dt_symb
x1, x2 = self.xx_symb
l, g, a, omega, gamma = self.pp_symb
# u0 = input force
u0 = self.uu_symb[0]
# create symbolic rhs functions
dx1_dt = x2
dx2_dt = -2 * gamma * x2 - (g / l + 1 / l * u0) * sp.sin(x1)
# put rhs functions into a vector
self.dxx_dt_symb = sp.Matrix([dx1_dt, dx2_dt])
return self.dxx_dt_symb
parameters.py
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 11 13:51:06 2021
@author: Jonathan Rockstroh
"""
import sys
import os
import numpy as np
import sympy as sp
import tabulate as tab
model_name = "Kapitzas_Pendulum"
# --------- CREATE SYMBOLIC PARAMETERS
pp_symb = [l, g, a, omega, gamma] = sp.symbols("l, g, a, omega, gamma", real=True)
# -------- CREATE AUXILIARY SYMBOLIC PARAMETERS
# (parameters, which shall not be numerical represented in the parameter tabular)
omega_0 = sp.Symbol("omega_0")
# --------- SYMBOLIC PARAMETER FUNCTIONS
# ------------ parameter values can be constant/fixed values OR
# ------------ set in relation to other parameters (for example: a = 2*b)
l_sf = 30 / 100
g_sf = 9.81
a_sf = 1 / 5 * l
omega_sf = 16 * omega_0
gamma_sf = 0.1 * omega_0
# List of symbolic parameter functions
pp_sf = [l_sf, g_sf, a_sf, omega_sf, gamma_sf]
# Set numerical values of auxiliary parameters
omega_0_nv = np.sqrt(g_sf / l_sf)
# List for Substitution
# -- Entries are tuples like: (independent symbolic parameter, numerical value)
pp_subs_list = [(l, l_sf), (omega_0, omega_0_nv)]
# OPTONAL: Dictionary which defines how certain variables shall be written
# in the tabular - key: Symbolic Variable, Value: LaTeX Representation/Code
# useful for example for complex variables: {Z: r"\underline{Z}"}
latex_names = {}
# ---------- CREATE BEGIN OF LATEX TABULAR
# Define tabular Header
# DON'T CHANGE FOLLOWING ENTRIES: "Symbol", "Value"
tabular_header = ["Parameter Name", "Symbol", "Value", "Unit"]
# Define column text alignments
col_alignment = ["left", "center", "left", "center"]
# Define Entries of all columns before the Symbol-Column
# --- Entries need to be latex code
col_1 = [
"Pendulum length",
"acceleration due to gravitation",
"Amplitude of Oscillation",
"Frequency of Oscillation",
"Dampening Factor",
]
# contains all lists of the columns before the "Symbol" Column
# --- Empty list, if there are no columns before the "Symbol" Column
start_columns_list = [col_1]
# Define Entries of the columns after the Value-Column
# --- Entries need to be latex code
col_4 = ["cm", r"$\frac{m}{s^2}$", "cm", "Hz", "Hz"]
# contains all lists of columns after the FIX ENTRIES
# --- Empty list, if there are no columns after the "Value" column
end_columns_list = [col_4]