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
from pyblocksim import *
# link to documentation with examples: https://ackrep-doc.readthedocs.io/en/latest/devdoc/contributing_data.html
def simulate():
"""
simulate the system model with blocks
:return: time, hysteresis output, input signal
"""
model = system_model.Model()
# input
uu = model.input_ramps
Result_block, hyst_in = model.get_blockfnc()
tt, states = blocksimulation(25, {hyst_in: uu}) # simulate
tt = tt.flatten()
bo = compute_block_ouptputs(states)
input_signal = [uu(t) for t in tt]
simulation_data = [tt, bo[Result_block], input_signal]
save_plot(simulation_data)
return simulation_data
def save_plot(simulation_data):
"""
plot your data and save the plot
:param simulation_data: simulation_data of system_model
:return: None
"""
# ---------start of edit section--------------------------------------
# plot of your data
fig1, axs = plt.subplots(nrows=2, ncols=1, figsize=(12.8, 9.6))
# print in axes top left
axs[0].plot(simulation_data[0], simulation_data[2], label="input")
axs[0].plot(simulation_data[0], simulation_data[1], label="hyst. output")
axs[0].set_xlabel("Time [s]") # x-label
axs[0].grid(1)
axs[0].legend()
# print in axes top right
axs[1].plot(simulation_data[2], simulation_data[1])
axs[1].set_xlabel("Input signal") # y-label
axs[1].set_ylabel("Hysteresis-output") # x-Label
axs[1].grid(1)
# ---------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 = [25, 2, 0]
# ---------end of edit section----------------------------------------
rc = ResultContainer(score=1.0)
simulated_final_state = [simulation_data[0][-1], simulation_data[1][-1], simulation_data[2][-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
import importlib
import sys, os
from pyblocksim import *
# 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 = 2
# ---------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
"""
u = sp.sin(4 * sp.pi * t)
return [u]
# ---------end of edit section----------------------------------------
return uu_rhs
def get_rhs_func(self):
msg = "This DAE model has no rhs func like ODE models."
raise NotImplementedError(msg)
def get_rhs_symbolic(self):
"""This model is not represented by the standard rhs equations."""
return False
def get_blockfnc(self):
"""
Calculate resulting blockfunction
:return: final block, input
"""
# ---------start of edit section--------------------------------------
s1, s2, y1, y2, T_storage = self.pp_symb # parameters
s1 = self.pp_dict[s1]
s2 = self.pp_dict[s2]
y1 = self.pp_dict[y1]
y2 = self.pp_dict[y2]
T_storage = self.pp_dict[T_storage]
_tanh_factor = 1e3
def step_factory(y0, y1, x_step):
"""
Factory to create continously approximated step functions
"""
# tanh maps R to (-1, 1)
# first map R to (0, 1)
# then map (0, 1) -> (y0, y1)
dy = y1 - y0
def fnc(x, module=sp):
return (module.tanh(_tanh_factor * (x - x_step)) + 1) / 2 * dy + y0
fnc.__doc__ = "approximated step function %f, %f, %f" % (y0, y1, x_step)
return fnc
# togehter these two are applied to the input
step1 = step_factory(-1, 0, s1)
step2 = step_factory(0, 1, s2)
# the third step-function limits the input of the PT1 between between 0 and 1
# the exact value for 63 % Percent
time_const_value = 1 - np.exp(-1)
step3 = step_factory(0, 1, time_const_value)
hyst_in, fb = inputs("hyst_in, fb") # overall input and internal feedback
# the sum of the two step-functions is basically a three-point-controller
SUM1 = Blockfnc(step1(hyst_in) + step2(hyst_in) + fb)
LIMITER = Blockfnc(step3(SUM1.Y))
PT1_storage = TFBlock(1 / (T_storage * s + 1), LIMITER.Y)
loop(PT1_storage.Y, fb)
# gain and offset for the output
Result_block = Blockfnc(PT1_storage.Y * (y2 - y1) + y1)
return [Result_block, hyst_in]
def input_ramps(self, t):
"""
input signal
:return: numerical value
"""
T1 = 10
T2 = 20
k1 = 1
k2 = 1
if t < 0:
return 0
elif t < T1:
return k1 * t
elif t < T2:
return k1 * T1 - k2 * (t - T1)
else:
return 0
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 = "Hysteresis system"
# ---------- create symbolic parameters
pp_symb = [s1, s2, y1, y2, T_storage] = sp.symbols("s1, s2, y1, y2, T_storage", 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)
s1_sf = 4
s2_sf = 8
y1_sf = 2
y2_sf = 11
T_storage_sf = 1e-4
# list of symbolic parameter functions
# tailing "_sf" stands for "symbolic parameter function"
pp_sf = [s1_sf, s2_sf, y1_sf, y2_sf, T_storage_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 = ["Symbol", "Value"]
# Define column text alignments
col_alignment = ["center", "left"]
# Define Entries of all columns before the Symbol-Column
# --- Entries need to be latex code
col_1 = []
# contains all lists of the columns before the "Symbol" Column
# --- Empty list, if there are no columns before the "Symbol" Column
start_columns_list = []
# Define Entries of the columns after the Value-Column
# --- Entries need to be latex code
col_4 = []
# contains all lists of columns after the FIX ENTRIES
# --- Empty list, if there are no columns after the "Value" column
end_columns_list = []