simulation.py
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
# link to documentation with examples: https://ackrep-doc.readthedocs.io/en/latest/devdoc/contributing_data.html
def simulate():
"""
simulate the system model with scipy.integrate.solve_ivp
:return: result of solve_ivp, might contains input function
"""
model = system_model.Model()
rhs_symb = model.get_rhs_symbolic()
print("Computational Equations:\n")
for i, eq in enumerate(rhs_symb):
print(f"dot_x{i+1} =", eq)
rhs = model.get_rhs_func()
# ---------start of edit section--------------------------------------
# initial state values
xx0 = [0, 0, 0, 0]
t_end = 5
tt = np.linspace(0, t_end, 10000)
simulation_data = solve_ivp(rhs, (0, t_end), xx0, t_eval=tt)
u = []
for i in range(len(simulation_data.t)):
u.append(model.uu_func(simulation_data.t[i], xx0)[0])
simulation_data.uu = u
# ---------end of edit section----------------------------------------
save_plot(simulation_data)
return simulation_data
def save_plot(simulation_data):
"""
plot your data and save the plot
access to data via: simulation_data.t array of time values
simulation_data.y array of data components
simulation_data.uu array of input values
:param simulation_data: simulation_data of system_model
:return: None
"""
# ---------start of edit section--------------------------------------
# plot of your data
fig1, axs = plt.subplots(nrows=5, ncols=1, figsize=(12.8, 9))
# print in axes top left
axs[0].plot(simulation_data.t, simulation_data.y[0])
axs[0].set_ylabel(r"$x_1 \, [m]$") # y-label
axs[0].grid()
axs[1].plot(simulation_data.t, simulation_data.y[1])
axs[1].set_ylabel(r"$x_2 \, [rad]$") # y-label
axs[1].grid()
axs[2].plot(simulation_data.t, simulation_data.y[2])
axs[2].set_ylabel(r"$x_3 \, [m/s]$") # y-label
axs[2].grid()
axs[3].plot(simulation_data.t, simulation_data.y[3])
axs[3].set_ylabel(r"$x_4 \, [rad/s]$") # y-label
axs[3].grid()
axs[4].plot(simulation_data.t, simulation_data.uu)
axs[4].set_ylabel(r"$u \, [Nm]$") # y-label
axs[4].set_xlabel("Time [s]") # x-label
axs[4].grid()
# ---------end of edit section----------------------------------------
plt.tight_layout()
save_plot_in_dir()
def evaluate_simulation(simulation_data):
"""
assert that the simulation results are as expected
:param simulation_data: simulation_data of system_model
:return:
"""
# ---------start of edit section--------------------------------------
# fill in final states of simulation to check your model
# simulation_data.y[i][-1]
expected_final_state = [50.46180780370284, -1.4107500985143662, 26.73084675897836, -0.044207917732603755]
# ---------end of edit section----------------------------------------
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
import sympy as sp
import symbtools as st
from symbtools import modeltools as mt
import importlib
import sys, os
# from ipydex import IPS, activate_ips_on_exception
from ackrep_core.system_model_management import GenericModel, import_parameters
# Import parameter_file
params = import_parameters()
# link to documentation with examples: https://ackrep-doc.readthedocs.io/en/latest/devdoc/contributing_data.html
class Model(GenericModel):
def initialize(self):
"""
this function is called by the constructor of GenericModel
:return: None
"""
# ---------start of edit section--------------------------------------
# Define number of inputs -- MODEL DEPENDENT
self.u_dim = 1
# Set "sys_dim" to constant value, if system dimension is constant
self.sys_dim = 4
# ---------end of edit section----------------------------------------
# check existence of params file
self.has_params = True
self.params = params
# ----------- SET DEFAULT INPUT FUNCTION ---------- #
# --------------- Only for non-autonomous Systems
def uu_default_func(self):
"""
define input function
:return:(function with 2 args - t, xx_nv) default input function
"""
# ---------start of edit section--------------------------------------
def uu_rhs(t, xx_nv):
"""
sequence of numerical input values
:param t:(scalar or vector) time
:param xx_nv:(vector or array of vectors) numeric state vector
:return:(list) numeric inputs
"""
if t < 1:
u_num_2 = 0
elif 1 < t < 1.5:
u_num_2 = -0.2
elif 1.5 < t < 3:
u_num_2 = 0.5
else:
u_num_2 = 0
return [u_num_2]
# ---------end of edit section----------------------------------------
return uu_rhs
# ----------- SYMBOLIC RHS FUNCTION ---------- #
def get_rhs_symbolic(self):
"""
define symbolic rhs function
:return: matrix of symbolic rhs-functions
"""
if self.dxx_dt_symb is not None:
return self.dxx_dt_symb
# ---------start of edit section--------------------------------------
x1, x2, x3, x4 = self.xx_symb # state components
m1, J1, J2, r, g = self.pp_symb # parameters
u1 = self.uu_symb[0] # inputs
q = sp.Matrix([[x1], [x2]])
xdot1, xdot2 = sp.symbols("xdot1, xdot2")
ex = sp.Matrix([1, 0])
ey = sp.Matrix([0, 1])
# coordinates of the centers of gravity
S1 = 0
S2 = x1 * mt.Rz(x2) * ex
# velocity of the ball
Sd2 = st.time_deriv(S2, q)
# kinetic energy
T_trans = (m1 * Sd2.T * Sd2) / 2 # translation energy
T_rot = (J1 * xdot2**2 + J2 * xdot2**2 + J2 * (xdot1 / r) ** 2) / 2 # rotation energy
T = T_trans[0] + T_rot
# potential Energie
V = m1 * g * S2[1]
external_forces = sp.Matrix([[0, u1]])
mod = mt.generate_symbolic_model(T, V, q, external_forces)
mod.calc_state_eq(simplify=False)
state_eq = mod.state_eq.subs([(xdot1, x3), (xdot2, x4)])
# ---------end of edit section----------------------------------------
return state_eq
parameters.py
import sys
import os
import numpy as np
import sympy as sp
import tabulate as tab
# link to documentation with examples: https://ackrep-doc.readthedocs.io/en/latest/devdoc/contributing_data.html
# set model name
model_name = "ball beam"
# ---------- create symbolic parameters
pp_symb = [m1, J1, J2, r, g] = sp.symbols("m1, J1, J2, r, g", real=True)
# ---------- create symbolic parameter functions
# parameter values can be constant/fixed values OR set in relation to other parameters (for example: a = 2*b)
m1_sf = 0.05
J1_sf = 0.02
J2_sf = 2e-6
r_sf = 0.01
g_sf = 10
# list of symbolic parameter functions
# tailing "_sf" stands for "symbolic parameter function"
pp_sf = [m1_sf, J1_sf, J2_sf, r_sf, g_sf]
# ---------- list for substitution
# -- entries are tuples like: (independent symbolic parameter, numerical value)
pp_subs_list = []
# OPTONAL: Dictionary which defines how certain variables shall be written
# in the table - key: Symbolic Variable, Value: LaTeX Representation/Code
# useful for example for complex variables: {Z: r"\underline{Z}"}
latex_names = {}
# ---------- Define LaTeX table
# Define table 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 = [
"mass of the ball",
"moment of inertia of the beam",
"moment of inertia of the ball",
"radius of the ball",
"acceleration due to gravity",
]
# 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 = ["kg", r"$kg \cdot m^2$", r"$kg \cdot m^2$", "m", r"$\frac{m}{s^2}$"]
# 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]