Casadi dae builder template and examples

Dockerfile

@@ -0,0 +1,48 @@
+# Stage 1: Build with nightly Rust
+FROM rustlang/rust:nightly AS builder
+
+WORKDIR /app
+
+COPY Cargo.toml Cargo.lock ./
+COPY src ./src
+COPY build.rs ./
+COPY modelica.par ./
+
+# Debug output
+RUN ls -R /app
+RUN cat Cargo.toml
+
+# Build with verbose logging
+RUN cargo build --release --verbose
+
+# Stage 2: Runtime image
+FROM python:3.12-slim
+
+
+# Install dependencies (for Rust build)
+RUN apt-get update && apt-get install -y \
+    build-essential \
+    cmake \
+    libssl-dev \
+    pkg-config \
+    bash \
+    && rm -rf /var/lib/apt/lists/*
+
+WORKDIR /tests
+
+COPY tests/notebooks/examples notebooks/examples
+COPY tests/notebooks/dependencies notebooks/dependencies
+COPY tests/models models
+COPY tests/templates templates
+COPY tests/instructions.md instructions.md
+
+RUN pip install --no-cache-dir -r notebooks/dependencies/requirements.txt
+RUN pip install notebooks/dependencies/casadi-3.7.0.dev+main-cp312-none-manylinux2014_x86_64.whl
+
+COPY --from=builder /app/target/release/rumoca /usr/local/bin/rumoca
+COPY --from=builder /app/modelica.par /usr/local/bin/modelica.par
+
+EXPOSE 8888
+
+ENTRYPOINT ["/bin/bash"]
+

tests/instructions.md

@@ -0,0 +1,25 @@
+Run Docker with port mapping (This allows jupyter notebook to run)
+sudo docker run -p 8888:8888 -it rumoca-container
+
+## Run Rumoca:
+Example 1: 
+
+rumoca --help
+
+Example 2 mass-spring-damper model to casadi:
+
+rumoca models/msd.mo -t templates/casadi_dae.jinja
+
+
+Example 3: save output to file:
+
+rumoca models/msd.mo -t templates/casadi_dae.jinja > msd_casadi_example.py
+
+
+
+## Run example notebooks to explore exported python models:
+jupyter lab --ip=0.0.0.0 --port=8888 --allow-root --NotebookApp.token=''
+
+Launch browser and go to: http://localhost:8888
+
+Explore jupyter notebooks.

tests/models/.ipynb_checkpoints/bouncing_ball-checkpoint.mo

@@ -0,0 +1,15 @@
+model BouncingBall "The 'classic' bouncing ball model"
+  parameter Real e=0.8 "Coefficient of restitution";
+  parameter Real h0=1.0 "Initial height";
+  Real h = 5.0 "Height";
+  Real v "Velocity";
+  Real z;
+equation
+  z = 2*h + v;
+  der(h) = v;
+  der(v) = -9.81;
+  when h<0 then
+    reinit(v, -e*pre(v));
+  end when;
+
+end BouncingBall;

tests/models/bouncing_ball.mo

@@ -1,14 +1,15 @@
 model BouncingBall "The 'classic' bouncing ball model"
   parameter Real e=0.8 "Coefficient of restitution";
   parameter Real h0=1.0 "Initial height";
-  Real h = 1.0 "Height";
+  Real h = 5.0 "Height";
   Real v "Velocity";
   Real z;
 equation
   z = 2*h + v;
-  v = der(h);
+  der(h) = v;
   der(v) = -9.81;
   when h<0 then
     reinit(v, -e*pre(v));
   end when;
+
 end BouncingBall;

tests/models/msd.mo

@@ -0,0 +1,12 @@
+model msd "Mass spring damper"
+  parameter Real k=1.0 "Coefficient of spring";
+  parameter Real c=1.0 "Coefficient of damper";
+  parameter Real m=1.0 "mass";
+  Real x = 0.0 "Position";
+  Real v "Velocity";
+  input Real u "Disturbance";
+equation
+  der(x) = v;
+  der(v) = -(k/m)*x - (c/m)*v - (1/m)*u;
+end msd;
+

tests/models/quadrotor.mo

@@ -10,18 +10,11 @@ model Quadrotor
     input Real e "elevator";
     input Real r "rudder";
     input Real t "throttle";
-    Real R_z "ground reaction force";
 equation
-    if h < 0 then
-        R_z = 10*h;
-    else
-        R_z = 0;
-    end if;
-
     // body forces
-    F_x = -(m*g - R_z)*sin(theta);
-    F_y = (m*g - R_z)*sin(phi)*cos(theta);
-    F_z = (m*g - R_z)*cos(phi)*cos(theta) -
+    F_x = -(m*g)*sin(theta);
+    F_y = (m*g)*sin(phi)*cos(theta);
+    F_z = (m*g)*cos(phi)*cos(theta) -
         (m_1.thrust + m_2.thrust + m_3.thrust + m_4.thrust);
 
     // body momments
@@ -68,9 +61,9 @@ equation
     der(psi) = (Q*sin(phi) + R*cos(phi))/cos(theta);
 
     // moment equations
-    Lambda*der(P) = J_xz*(J_x - J_y + J_z)*P*Q - (J_z*(J_z - J_y) + J_xz*J_xz)*Q*R + J_z*M_x + J_xz*M_z;
-    J_y*der(Q) = (J_z - J_x)*P*R - J_xz*(P*P - R*R) + M_y;
-    Lambda*der(R) = ((J_x - J_y)*J_x + J_xz*J_xz)*P*Q - J_xz*(J_x - J_y + J_z)*Q*R + J_xz*M_x + J_x*M_z;
+    der(P) = (J_xz*(J_x - J_y + J_z)*P*Q - (J_z*(J_z - J_y) + J_xz*J_xz)*Q*R + J_z*M_x + J_xz*M_z)/Lambda;
+    der(Q) = ((J_z - J_x)*P*R - J_xz*(P*P - R*R) + M_y)/J_y;
+    der(R) = (((J_x - J_y)*J_x + J_xz*J_xz)*P*Q - J_xz*(J_x - J_y + J_z)*Q*R + J_xz*M_x + J_x*M_z)/Lambda;
 
 
 end RigidBody6DOF;

tests/notebooks/bouncing_ball_ca.py

@@ -0,0 +1,134 @@
+"""
+Automatically generated by Rumoca
+"""
+import casadi as ca
+import numpy as np
+sin = ca.sin
+cos = ca.cos
+tan = ca.tan
+
+class Model:
+    """
+    Flattened Modelica Model
+    """
+
+    def __init__(self, model_name):
+        # ============================================
+        # Initialize
+
+
+        # Create new dae
+        dae = ca.DaeBuilder(model_name)
+        # ============================================
+        # Declare time
+        time = dae.add('time', 'independent')
+        # If equations update time
+        dt = 0.01
+
+        # ============================================
+        # Declare u
+
+        # ============================================
+        # Declare p
+
+        e = dae.add('e', 'parameter', 'tunable', dict(start = 0.8))
+        h0 = dae.add('h0', 'parameter', 'tunable', dict(start = 1.0))
+        # ============================================
+        # Declare c# ============================================
+        # Declare cp
+
+        # ============================================
+        # Declare x
+
+        h = dae.add('h', dict(start = 5.0))
+        v = dae.add('v', dict(start = 0.0))
+        # ============================================
+        # Declare m
+
+        # ============================================
+        # Declare y
+
+        z = dae.add('z', dict(start = 0.0))
+        # ============================================
+        # Declare z
+
+
+
+
+        # ============================================
+        # Declare pre_x
+        pre_h = dae.pre(h)
+        pre_v = dae.pre(v)
+        # ============================================
+        # Declare pre_m
+        # ============================================
+        # Declare pre_z
+        # ============================================
+        # Define Condition Update Function: fc
+
+        def c0():
+            return (h < 0.0)
+
+        # ============================================
+        # Define reset functions: fr
+        def c0_fr():
+            return dae.reinit("v",-((e * pre_v)))
+        dae.when(c0(), [c0_fr()])
+        # ============================================
+        # Declare x_dot
+        der_h = dae.der(h)
+        der_v = dae.der(v)
+        # ============================================
+        def if_else_builder(s, builder, terminal_state):
+            state = terminal_state
+            for i, (cond, value) in enumerate(reversed(builder)):
+                state = ca.if_else(s == int(cond[1:]), value, state)
+            return state
+
+        def if_else_builder2(builder, terminal_state):
+            state = terminal_state
+            for cond, value in reversed(builder):
+                state = ca.if_else(cond, value, state)
+            return state
+        # ============================================
+
+        def add_expression(dictionary, var, expression):
+            if var not in dictionary.keys():
+                dictionary[var] = []
+            dictionary[var].append(expression)
+
+            return dictionary
+        # Define Continous Update Function: fx
+        dae.eq(z, ((2.0 * h) + v))
+        dae.eq(der_h, v)
+        dae.eq(der_v, -(9.81))
+
+
+
+        dae.sort('w')
+        self.dae = dae
+    def display(self):
+        self.dae.disp(True)
+
+    def simulate(self, t0, tf, dt, x0=None, p0=None, f_u=None, max_events=100):
+        """
+        Simulate the modelica model
+        """
+        if p0 is None:
+            p0 = self.dae.start(self.dae.p())
+
+        if x0 is None:
+            x0 = self.dae.start(self.dae.x())
+
+        tgrid = np.arange(t0, tf, dt)
+        simopts = dict(transition = self.dae.transition(), verbose = False,
+            event_tol = 1e-12, max_events = 100000, max_event_iter = 2000)
+
+        sim = ca.integrator('sim', 'cvodes', self.dae.create(), 0, tgrid, simopts)
+
+        if f_u is None:
+            simres = sim(x0 = x0, p = p0)
+        else:
+            simres = sim(x0 = x0, p = p0, u=f_u)
+
+        return tgrid, simres

tests/notebooks/dependencies/requirements.txt

@@ -0,0 +1,7 @@
+control==0.10.1
+matplotlib==3.10.1
+pandas==2.2.3
+scipy==1.15.2
+sympy==1.14.0
+jupyterlab==4.0.9
+

tests/notebooks/examples/.ipynb_checkpoints/bouncing_ball_casadi-checkpoint.ipynb

@@ -0,0 +1,87 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The autoreload extension is already loaded. To reload it, use:\n",
+      "  %reload_ext autoreload\n",
+      "lexer invalid token InvalidToken\n"
+     ]
+    }
+   ],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "! rumoca ../../models/bouncing_ball.mo -t ../../templates/casadi_dae.jinja > bouncing_ball_ca.py"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "SyntaxError",
+     "evalue": "invalid syntax (bouncing_ball_ca.py, line 1)",
+     "output_type": "error",
+     "traceback": [
+      "Traceback \u001b[0;36m(most recent call last)\u001b[0m:\n",
+      "\u001b[0m  File \u001b[1;32m~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3579\u001b[0m in \u001b[1;35mrun_code\u001b[0m\n    exec(code_obj, self.user_global_ns, self.user_ns)\u001b[0m\n",
+      "\u001b[0;36m  Cell \u001b[0;32mIn[7], line 1\u001b[0;36m\n\u001b[0;31m    import bouncing_ball_ca\u001b[0;36m\n",
+      "\u001b[0;36m  File \u001b[0;32m~/Research/devlopment/rumoca/tests/notebooks/examples/bouncing_ball_ca.py:1\u001b[0;36m\u001b[0m\n\u001b[0;31m    lexer invalid token InvalidToken\u001b[0m\n\u001b[0m          ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
+     ]
+    }
+   ],
+   "source": [
+    "import bouncing_ball_ca\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "model = bouncing_ball_ca.Model(\"bouncing_ball\")\n",
+    "tgrid, res = model.simulate(t0=0, tf=8, dt=0.01)\n",
+    "\n",
+    "plt.plot(tgrid, res['xf'].T, label=model.dae.x())\n",
+    "\n",
+    "# Set axis label and tick font sizes\n",
+    "plt.xlabel(\"Time (s)\", fontsize=16)\n",
+    "plt.ylabel(\"State\", fontsize=16)\n",
+    "plt.tick_params(axis='both', which='major', labelsize=14)\n",
+    "\n",
+    "# Make legend font larger\n",
+    "plt.legend(fontsize=14)\n",
+    "plt.grid()\n",
+    "\n",
+    "# Save without cutting off labels\n",
+    "plt.savefig(\"casadi_bouncing_ball.pdf\", bbox_inches='tight')\n",
+    "\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}