Changes related to using libcosimpy with the ECCO algorithm. Test is included in test_oscillator_fmu

Author: eisDNV

examples/ForcedOscillator.xml

@@ -7,9 +7,13 @@
       <Simulator name="drv" source="DrivingForce.fmu" stepSize="0.01" />
    </Simulators>
    <Connections>
-      <VariableConnection>
-         <Variable simulator="drv" name="f[2]" />
-         <Variable simulator="osc" name="f[2]" />
+      <VariableConnection powerBond="force-velocity">
+         <Variable simulator="drv" name="f[2]" causality="output"/>
+         <Variable simulator="osc" name="f[2]" causality="input"/>
+      </VariableConnection>
+      <VariableConnection powerBond="force-velocity">
+         <Variable simulator="osc" name="v[2]" causality="output"/>
+         <Variable simulator="drv" name="v_osc[2]" causality="input"/>
       </VariableConnection>
    </Connections>
    <EccoConfiguration>

examples/bouncing_ball_3d.py

@@ -130,7 +130,7 @@ def _interface(self, name: str, start: str | float | tuple) -> Variable:
                 variability="continuous",
                 initial="exact",
                 start=start,
-                rng=((0, "100 m"), None, (0, "10 m")),
+                rng=((0, "100 m"), None, (0, "39.371 inch")),
             )
         elif name == "speed":
             return Variable(

examples/driving_force_fmu.py

@@ -21,20 +21,21 @@ def func(time: float, ampl: float = 1.0, omega: float = 0.1, d_omega: float = 0.
     """
     if d_omega == 0.0:
         return np.array((0, 0, ampl * sin(omega * time)), float)
-    else:           
-        return np.array((0, 0, ampl * sin((omega + d_omega*time) * time)), float)
+    else:
+        return np.array((0, 0, ampl * sin((omega + d_omega * time) * time)), float)
 
 
 class DrivingForce(Model):
     """A driving force in 3 dimensions which produces an ouput per time and can be connected to the oscillator.
 
+    Note1: the FMU model is made directly (without a basic python class model), which is not recommended!
+    Note2: the speed of the connected oscillator is added as additional connector.
+    Since the driving is forced, the input speed is ignored, but it is needed for the ECCO algorithm (power bonds).
+
     Note: This completely replaces DrivingForce (do_step and other functions are not re-used).
 
     Args:
         func (callable)=func: The driving force function f(t).
-            Note: The func can currently not really be handled as parameter and must be hard-coded here (see above).
-            Soon to come: Model.build() function which honors parameters, such that function can be supplied from
-            outside and the FMU can be re-build without changing the class.
     """
 
     def __init__(
@@ -51,7 +52,13 @@ def __init__(
             "Siegfried Eisinger",
             **kwargs,
         )
-        # interface Variables
+        # interface Variables. We define first their values, to help pyright, since the basic model is missing
+        self.ampl: float = ampl
+        self.freq: float = freq
+        self.d_freq: float = d_freq
+        self.function: Callable = function
+        self.f = np.array((0, 0, 0), float)
+        self.v_osc = np.array((0, 0, 0), float)
         self._ampl = Variable(self, "ampl", "The amplitude of the force in N", start=ampl)
         self._freq = Variable(self, "freq", "The frequency of the force in 1/s", start=freq)
         self._d_freq = Variable(self, "d_freq", "Change of frequency of the force in 1/s**2", start=d_freq)
@@ -64,15 +71,25 @@ def __init__(
             variability="continuous",
             start=np.array((0, 0, 0), float),
         )
+        self._v_osc = Variable(
+            self,
+            "v_osc",
+            "Input connector for the speed of the connected element in m/s",
+            causality="input",
+            variability="continuous",
+            start=np.array((0, 0, 0), float),
+        )
 
     def do_step(self, current_time: float, step_size: float):
-        self.f = self.func(current_time+step_size)
+        self.f = self.func(current_time + step_size)
         return True  # very important!
 
     def exit_initialization_mode(self):
         """Set internal state after initial variables are set."""
-        self.func = partial(self.function,
-                            ampl=self.ampl,
-                            omega=2 * pi * self.freq,
-                            d_omega= 0.0 if self.d_freq == 0.0 else 2 * pi * self.d_freq)
+        self.func = partial(
+            self.function,
+            ampl=self.ampl,
+            omega=2 * pi * self.freq,
+            d_omega=0.0 if self.d_freq == 0.0 else 2 * pi * self.d_freq,
+        )
         logger.info(f"Initial settings: ampl={self.ampl}, freq={self.freq}, d_freq={self.d_freq}")

examples/oscillator.py

@@ -33,7 +33,7 @@ def __init__(
         c: tuple[float, float, float] | tuple[str, str, str] = (0.0, 0.0, 0.0),
         m: float = 1.0,
         tolerance: float = 1e-5,
-        f_func: Callable|None = None
+        f_func: Callable | None = None,
     ):
         self.k = np.array(k, float)
         self.c = np.array(c, float)
@@ -51,8 +51,8 @@ def ode_func(self, t: float, y: np.ndarray, i: int, f: float) -> np.ndarray:
         if self.f_func is None:
             if f != 0:
                 res += np.array((f, 0), float)
-        elif i==2: # only implemented for z
-            res += np.array((self.f_func(t)[i], 0), float)            
+        elif i == 2:  # only implemented for z
+            res += np.array((self.f_func(t)[i], 0), float)
         return res
 
     def do_step(self, current_time: float, step_size: int | float) -> bool:

src/component_model/variable.py

@@ -68,7 +68,7 @@ def __str__(self):
         return txt
 
     def parse_quantity(self, quantity: PyType, ureg: UnitRegistry, typ: type | None = None) -> PyType:
-        """Disect the provided quantity in terms of magnitude and unit, if provided as string.
+        """Parse the provided quantity in terms of magnitude and unit, if provided as string.
         If another type is provided, dimensionless units are assumed.
 
         Args:

tests/test_axle_fmu.py

@@ -146,7 +146,7 @@ def test_use_fmu(axle_fmu: Path, show: bool):
 
 
 if __name__ == "__main__":
-    retcode = 0#pytest.main(["-rP -s -v", __file__])
+    retcode = 0  # pytest.main(["-rP -s -v", __file__])
     assert retcode == 0, f"Return code {retcode}"
     import os
 

tests/test_bouncing_ball_3d.py

@@ -420,6 +420,8 @@ def get_status(sim):
     assert values[6] == 1.5, "Initial setting did not work"
     assert values[5] == -0.015, "Initial setting did not have the expected effect on speed"
 
+    manipulator.slave_real_values(ibb, [0, 1, 2], [1.0, 2.0, 3.0])
+
 
 #     values = observer.last_real_values(0, list(range(11)))
 #     print("VALUES2", values)
@@ -428,6 +430,7 @@ def get_status(sim):
 
 #    sim.simulate_until(target_time=3e9)
 
+
 def test_from_fmu(bouncing_ball_fmu):
     assert bouncing_ball_fmu.exists(), "FMU not found"
     model = model_from_fmu(bouncing_ball_fmu)
@@ -450,13 +453,13 @@ def test_from_fmu(bouncing_ball_fmu):
 
 
 if __name__ == "__main__":
-    retcode = 0  # pytest.main(["-rA", "-v", "--rootdir", "../", "--show", "False", __file__])
+    retcode = pytest.main(["-rA", "-v", "--rootdir", "../", "--show", "False", __file__])
     assert retcode == 0, f"Non-zero return code {retcode}"
     import os
 
     os.chdir(Path(__file__).parent / "test_working_directory")
     # test_bouncing_ball_class(show=False)
-    test_make_bouncing_ball(_bouncing_ball_fmu())
+    # test_make_bouncing_ball(_bouncing_ball_fmu())
     # test_use_fmu(_bouncing_ball_fmu(), True)
     # test_from_fmu( _bouncing_ball_fmu())
-    # test_from_osp( _bouncing_ball_fmu())
+    # test_from_osp(_bouncing_ball_fmu())

tests/test_oscillator.py

@@ -1,9 +1,9 @@
 from functools import partial
 from math import atan2, cos, exp, pi, sin, sqrt
+from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
-from pathlib import Path
 
 
 def arrays_equal(res: tuple[float, ...] | list[float], expected: tuple[float, ...] | list[float], eps=1e-7):
@@ -29,9 +29,9 @@ def do_show(time: list, z: list, v: list, compare1: list | None = None, compare2
 
 def force(t: float, ampl: float = 1.0, omega: float = 0.1, d_omega: float = 0.0):
     if d_omega == 0.0:
-        return np.array((0, 0, ampl * sin(omega * t)), float) # fixed frequency
+        return np.array((0, 0, ampl * sin(omega * t)), float)  # fixed frequency
     else:
-        return np.array((0, 0, ampl * sin((omega + d_omega*t) * t)), float) # frequency sweep
+        return np.array((0, 0, ampl * sin((omega + d_omega * t) * t)), float)  # frequency sweep
 
 
 def forced_oscillator(
@@ -118,6 +118,7 @@ def run_oscillation_z(
 
     return (osc, times, z, v)
 
+
 def sweep_oscillation_z(
     k: float,
     c: float,
@@ -136,13 +137,9 @@ def sweep_oscillation_z(
     and return the oscillator object and the time series for z-position and z-velocity."""
 
     from examples.oscillator import Oscillator
-    f_func = f_func=partial(force, ampl=ampl, omega=0.0, d_omega=d_omega)
-    osc = Oscillator(k=(1.0, 1.0, k),
-                     c=(0.0, 0.0, c),
-                     m=m,
-                     tolerance=tol,
-                     f_func = f_func
-                     )
+
+    f_func = f_func = partial(force, ampl=ampl, omega=0.0, d_omega=d_omega)
+    osc = Oscillator(k=(1.0, 1.0, k), c=(0.0, 0.0, c), m=m, tolerance=tol, f_func=f_func)
     osc.x[2] = x0  # set initial z value
     osc.v[2] = v0  # set initial z-speed
     times, z, v, f = [], [], [], []
@@ -264,70 +261,75 @@ def area(x: list[float], y: list[float]):
     assert abs(y[-1] - 4.0) < 1e-12, f"Found {y[-1]}"
     osc, x, y = run_2d(x0=(1.0, 0.0, 0.0), v0=(0.0, 1.0, 0.0), k=(1.0, 1.0 / 15.8, 0), end=20 * np.pi)
     if show:
-        show_2d(x,y)
+        show_2d(x, y)
+
 
 def test_sweep_oscillator(show: bool = True):
     """A forced oscillator where the force frequency is changed linearly as d_omega*time.
     The test demonstrates that a monolithic simulation provides accurate results in all ranges of the force frequency.
     Co-simulating the oscillator and the force, this does not work.
     """
-    osc, times0, z0, v0, f0 = sweep_oscillation_z( k=1.0, 
-                                               c=0.1,
-                                               m=1.0,
-                                               ampl=1.0,
-                                               d_omega=0.1,
-                                               x0=0.0,
-                                               v0=0.0,
-                                               dt=0.1, # 'ground truth', small dt
-                                               end=100.0,
-                                               tol=1e-3
-                                               )
-    with open(Path.cwd() / "oscillator_sweep0.dat", 'w') as fp:
+    osc, times0, z0, v0, f0 = sweep_oscillation_z(
+        k=1.0,
+        c=0.1,
+        m=1.0,
+        ampl=1.0,
+        d_omega=0.1,
+        x0=0.0,
+        v0=0.0,
+        dt=0.1,  # 'ground truth', small dt
+        end=100.0,
+        tol=1e-3,
+    )
+    with open(Path.cwd() / "oscillator_sweep0.dat", "w") as fp:
         for i in range(len(times0)):
-            fp.write( f"{times0[i]}\t{z0[i]}\t{v0[i]}\t{f0[i]}\n")
-            
+            fp.write(f"{times0[i]}\t{z0[i]}\t{v0[i]}\t{f0[i]}\n")
+
     if show:
-        freq = [0.1*t/2/np.pi for t in times0]
+        freq = [0.1 * t / 2 / np.pi for t in times0]
         fig, ax = plt.subplots()
         ax.plot(freq, z0, label="z0(t)")
         ax.plot(freq, v0, label="v0(t)")
         ax.plot(freq, f0, label="F0(t)")
         ax.legend()
         plt.show()
-        
-    osc, times, z, v, f = sweep_oscillation_z( k=1.0,
-                                               c=0.1,
-                                               m=1.0,
-                                               ampl=1.0,
-                                               d_omega=0.1,
-                                               x0=0.0,
-                                               v0=0.0,
-                                               dt=1, # dt similar to resonance frequency
-                                               end=100.0,
-                                               tol=1e-3
-                                               )
+
+    osc, times, z, v, f = sweep_oscillation_z(
+        k=1.0,
+        c=0.1,
+        m=1.0,
+        ampl=1.0,
+        d_omega=0.1,
+        x0=0.0,
+        v0=0.0,
+        dt=1,  # dt similar to resonance frequency
+        end=100.0,
+        tol=1e-3,
+    )
     i0 = 0
-    for i in range(len(times)): # demonstrate that the results are accurate, even if dt is large
+    for i in range(len(times)):  # demonstrate that the results are accurate, even if dt is large
         t = times[i]
-        while abs(times0[i0]-t) > 1e-10:
+        while abs(times0[i0] - t) > 1e-10:
             i0 += 1
             assert times0[i0] - t < 0.1, f"Time entry for time {t} not found in times0"
-            
+
         assert abs(z0[i0] - z[i]) < 2e-2, f"Time {t}. Found {z0[i0]} != {z[i]}"
         assert abs(v0[i0] - v[i]) < 2e-2, f"Time {t}. Found {v0[i0]} != {v[i]}"
-        
+
     if show:
         fig, ax = plt.subplots()
         ax.plot(times0, z0, label="z0(t)")
         ax.plot(times, z, label="z(t)")
         ax.legend()
         plt.show()
 
+
 if __name__ == "__main__":
-    retcode = 0  # pytest.main(["-rA", "-v", "--rootdir", "../", "--show", "False", __file__])
+    retcode = pytest.main(["-rA", "-v", "--rootdir", "../", "--show", "False", __file__])
     assert retcode == 0, f"Non-zero return code {retcode}"
     import os
+
     os.chdir(Path(__file__).parent.absolute() / "test_working_directory")
     # test_oscillator_class(show=True)
     # test_2d(show=True)
-    test_sweep_oscillator()
+    # test_sweep_oscillator()