DAE Examples and Last Point Strategy

README.md

@@ -48,7 +48,7 @@ We set ODE solver and its options as usual. If not set, default values are used.
       initIFDIFF();  % Initialise the paths for ifdiff -- needed only once per Matlab session
       integrator = @ode45; 
       odeoptions = odeset('AbsTol', 1e-14, 'RelTol', 1e-12);
-      datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', func2str(integrator), 'options', odeoptions);
+      datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', integrator, 'options', odeoptions);
    ```
 
 

toolbox/examples/daeExample/daeExampleRHS.m

@@ -0,0 +1,13 @@
+function f = daeExampleRHS(~, x, p)
+
+    f = zeros(2,1);
+    % algebraic constraint
+    z = x(1) + x(2);
+    f(2) = z;
+    % differential variables
+    if x(2) < p
+        f(1) = x(2);
+    else
+        f(1) = 0;
+    end
+end

toolbox/examples/daeExample/daeExample_README.md

@@ -0,0 +1,112 @@
+<script
+  src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
+  type="text/javascript">
+</script>
+
+
+# Switched DAE Example
+
+## Introduction
+
+
+A **differential algebraic equation (DAE)** is a system of equations involving an unknown function and its derivatives together with algebraic equations, which impose constraints on the variables. Geometrically, the algebraic equations define a constraint manifold in the state space, and the solution trajectory of any DAE initial-value problem must remain on this manifold.
+
+In applications, DAEs can be combined with state-dependent events (e.g. contact/no-contact, switching circuits), which leads to systems whose dynamics change when certain state conditions are met. IFDIFF can solve such switched DAEs.
+
+Let's take a look at the following example for $t \in [0,n]$ where  $n \in \mathbb{N}$:
+
+$$(D) \quad \begin{cases} \dot{x}_ 1 = f_1(x_2) = \begin{cases} x_2, \quad \text{if}  \quad x_2 < p \\ 0 \quad \text{if} \quad x_2 \geq p \end{cases}  \\ f_2(x) = x_1 + x_2 = 0  \\ x_0 = (1,-1)^T, p \in \mathbb{R} \end{cases}$$
+
+
+## Analytical solution
+
+First we derive the analytical solution of $(D)$. We see that the DAE is of index 1 since the algebraic constraint can be differentiated once to achieve $x_1 + x_2 = 0 \Rightarrow \dot{x_1} + \dot{x_2} = 0 \iff \dot{x_2} = - \dot{x_1}$ which then gives us the ODE system:
+
+$$(D_{\text{ODE}})\begin{cases} \dot{x}_ 1  = f_1(x_2) = \begin{cases} x_2, \quad \text{if} \hspace{0.2cm} x_2 < p\\ 0 \hspace{0.2cm} \text{if} \quad x_2 \geq p \end{cases}  \\ f_2(x) = - \dot{x_1} = \begin{cases} -x_2 \quad \text{if} \quad x_2 < p\\ 0 \quad \text{if} \quad x_2 \geq p \end{cases}   \\ x_0 = (1,-1)^T, p \in \mathbb{R} \end{cases}$$
+
+The system $(D_{\text{ODE}})$ will only exhibit switching behaviour at $p \in (-1,0)$ since:
+
+1. Let $ p \leq -1 $ : 
+At $ t = 0 $ we have $x_2(0)= -1 \geq p$. Therefore $\dot{x}_ 2(t)=0$, meaning we stay in the second branch of $f_1$ and no switch occurs.
+
+2. Let $p \geq 0 $ : 
+At $ t = 0 $ we have $x_2(0)= -1 < p$.Therefore we only have $\dot{x_2} = - \dot{x_1} \Rightarrow x_1(t) = e^{-t}$ and $x_2(t) = -e^{-t} \hspace{0.2cm} \forall t \in I$. Since $-e^{-t} < 0 \leq p$, not switch occurs in this situation either.
+
+3. So, let $p \in (-1, 0)$:
+In this situation the system starts in case 2, so exponential growth in the 1st and exponential decay in the 2nd component.
+We now determine the switching point $t_s \in (0,2)$. The switching point satisfies $x_2(t_s)=p$, so $x_2(t_s)=p \iff -e^{-t_s} = p \iff t_s = -\text{ln}(-p) $. (Remember that $-p$ is positive!)
+
+We also notice from the formula $t_s = -\text{ln}(-p)$ that the switching point $t_s$ will be larger as $p$ gets smaller. This means, depending on the parameter, we need to choose a fitting time horizon $[0,n]$.
+
+## Solution with IFDIFF
+
+In practice, DAE systems are not solved by converting them to an ODE systems. Instead, MATLAB offers the solver `ode15s` for solving DAEs which we will use IFDIFF with.
+
+### Step 1: Right Hand Side
+We code the Right Hand Side (RHS) as follows.
+
+```
+function f = daeExampleRHS(~, x, p)
+f = zeros(2,1);
+
+% algebraic constraint
+z = x(1) + x(2);
+f(2) = z;
+
+% differential variables
+if x(2) < p
+    f(1) = x(2);
+else
+    f(1) = 0;
+end
+```
+
+### Step 2: Setup & Integration
+For the main script, we set up a consistent initial value, e.g. $x_0 = (1,-1)^T$ and a mass matrix $M$ such that the algebraic constraint is set to 0 while the differential variables remain as coded in the RHS.
+We choose a parameter $ p \in (-1,0) $, e.g. $p = -0.2$ and a suitable time horizon, e.g. $[0 5]$.
+
+```
+integrator = @ode15s;
+x0         = [1; -1];
+tspan      = [0 5];
+M          = [1 0; 0 0];
+p          = -0.2;
+
+opts_ifdiff = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-6);
+opts_plain  = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9, 'RelTol', 1e-6);
+
+datahandle = prepareDatahandleForIntegration('daeExampleRHS', 'integrator', integrator, 'options', opts_ifdiff);
+sol_ifdiff = solveODE(datahandle, tspan, x0, p);
+sol_plain  = integrator(@(t, x) daeExampleRHS(t, x, p), tspan, x0, opts_plain);
+```
+
+### Step 3: Visualising
+
+To look at our solution and compare them to the plain `ode15s` we can plot both in one plot as follows.
+
+```
+fig1 = figure(01);
+hold on
+IFDIFF_plot_1 = plot(sol_ifdiff.x, sol_ifdiff.y, 'ro--', 'DisplayName', 'IFDIFF'); 
+Plain_plot_1  = plot(sol_plain.x, sol_plain.y, 'ko-', 'DisplayName', 'plain ode15s');
+Switch_plot   = xline(sol_ifdiff.switches, 'g', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+legend([Plain_plot_1(1), IFDIFF_plot_1(1), Switch_plot]);
+hold off
+```
+
+
+
+
+![](plots_daeExample/plot1.png)
+
+We notice that the integrator strategy results in small steps here for the plain solver as well as IFDIFF. This is because ode15 is a multi-step method and can thus not be changed.
+However, if we take a closer look, we see, that the integration with IFDIFF is accurate around the switching point. 
+
+![](plots_daeExample/plot1_close2.png)
+
+
+
+## Additional Content
+
+To further investigate this example, take a look at the files `daeExample_main.m` and `daeExampleRHs.m`. 
+Additionally, go to the folder `rlcExample` to learn about solving DAEs with IFDIFF in a modelling example.

toolbox/examples/daeExample/daeExample_main.m

@@ -0,0 +1,26 @@
+% DAE Example Main
+% for more details, see daeExample_README.md
+
+%% Setup and integration
+integrator = @ode15s;
+x0         = [1; -1];
+tspan      = [0 5];
+M          = [1 0; 0 0];
+p          = -0.2;
+
+opts_ifdiff = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-6);
+opts_plain  = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9, 'RelTol', 1e-6);
+
+datahandle = prepareDatahandleForIntegration('daeExampleRHS', 'integrator', integrator, 'options', opts_ifdiff);
+sol_ifdiff = solveODE(datahandle, tspan, x0, p);
+sol_plain  = integrator(@(t, x) daeExampleRHS(t, x, p), tspan, x0, opts_plain);
+
+%% Plots
+clf;
+fig1 = figure(01);
+hold on;
+IFDIFF_plot_1 = plot(sol_ifdiff.x, sol_ifdiff.y, 'ro--', 'DisplayName', 'IFDIFF'); 
+Plain_plot_1  = plot(sol_plain.x, sol_plain.y, 'ko-', 'DisplayName', 'plain ode15s');
+Switch_plot   = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+legend([Plain_plot_1(1), IFDIFF_plot_1(1), Switch_plot]);
+hold off;

toolbox/examples/rlcExample/rlcExample_README.md

@@ -0,0 +1,119 @@
+<script
+  src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
+  type="text/javascript">
+</script>
+
+# RLC Circuit (Voltage-controlled switch)
+
+This example models an **RLC electrical circuit** with a fuse using a differential-algebraic equation (DAE) formulation. 
+The system exhibits state-dependent switching between two modes depending on the capacitor voltage $V_C$ relative to a threshold $V_{\text{th}}$.
+
+## Physical Description
+
+The circuit consists of:
+
+- an inductor $L$
+
+- a capacitor $C$
+
+- two possible resistanive elements $R_1$ and $R_2$
+
+- a DC supply voltage $V_s$
+
+When the capacitor voltage exceeds a threshold $V_{\text{th}}$, the circuit switches from high to low resistance.
+Such a behavior may be observed​
+
+## Model
+
+The system states are $ x = (i_L, V_C, i_C)^T$. The inductor $L$ follows
+Kirchhoff's Voltage Law for the supply $V_s=Ri_L+L\dot{i}_L + V_C $ ($R$ is $R_1$ or $R_2$ depending on the case).
+The capacitator $C$ follows the capacitor-current-voltage relation $i_C=C\dot{V}_C$ where $V_C$ is the potential difference between
+the capacitor's positive and negative plates.
+
+The algebraic constraint in out model is Kirchhoff's Current Law at the inductor-capacitor node:
+Inductor current equals capacitor current since they are in series in the circuit.
+
+With this, the DAE system is given by:
+
+$ (\text D_{\text{RLC}} ) \quad \begin{cases} \dot{i}_L = \begin{cases} \frac{V_s - R_1 i_L - V_C}{L} \quad \text{if}  \quad V_C > V_{\text{th}} \\ \frac{V_s - R_2 i_L - V_C}{L} \quad \text{if} \quad V_C \leq V_{\text{th}} \end{cases} \\ \dot{V}_C = \frac{i_C}{C} \\ i_L−i_C = 0 \end{cases}$
+
+The system is an index 1 DAE since differentiating yields the second-order ODE $L\ddot{i}_L+R\dot{i}L+\frac{1}{C}i_L=0$.
+
+
+## Solution with IFDIFF
+
+The parameters are stored in a vector `p = [L; R1; R2; C; Vs; Vth]` since we need to keep the input structure `(t,x,p)` for IFDIFF.
+
+```
+function dx = rlcRHS(~,x,p)
+    dx  = zeros(3,1);
+    L   = p(1); R1  = p(2); R2  = p(3);
+    C   = p(4); Vs  = p(5); Vth = p(6);
+    iL  = x(1); VC  = x(2); iC  = x(3);
+
+    % Mode 1: High resistance (fuse intact)
+    if VC > Vth
+        dx(1) = (Vs - R1*iL - VC)/L;
+    else
+    % Mode 2: Low resistance (fuse blown)
+        dx(1) = (Vs - R2*iL - VC)/L;
+    end 
+    dx(2) = iC/C;
+    dx(3) = iL - iC; % algebraic constraint
+end
+```
+
+Now, we need to set up a consistent initial value, e.g. $x_0 =(0,0,0)^T$, a suitable time span and a mass matrix $M$ such that the algebraic constraint is set to 0 while the differential variables remain as coded in the RHS above.
+
+```
+integrator = @ode15s;
+M = diag([1,1,0]);
+x0 = [0; 0; 0];
+tspan = [0 10];
+p = [L; R1; R2; C; Vs; Vth];
+
+opts_ifdiff = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5);
+opts_plain  = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5);
+
+datahandle = prepareDatahandleForIntegration('rlcRHS', 'integrator', integrator, 'options', opts_ifdiff);
+
+sol_ifdiff = solveODE(datahandle, tspan, x0, p);
+sol_plain  = integrator(@(t, x) rlcRHS(t, x, p), tspan, x0, opts_plain);
+
+```
+To compare both solutions, we can plot them and mark the switch.
+
+```
+clf;
+fig1 = figure(01);
+
+subplot(2,1,1);
+hold on;
+Plot_ifdiff_1   = plot(sol_ifdiff.x, sol_ifdiff.y(1,:), 'ro--', 'DisplayName', 'IFDIFF'); 
+Plot_plain_1    = plot(sol_plain.x, sol_plain.y(1,:), 'k.-', 'DisplayName', 'plain ode15s');
+Switch_plot_1   = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+ylabel('i_L (A)');
+xlabel('Time (s)');
+legend();
+hold off;
+
+subplot(2,1,2);
+hold on
+Plot_ifdiff_2  = plot(sol_ifdiff.x, sol_ifdiff.y(2,:), 'ro--', 'DisplayName', 'IFDIFF' );
+Plot_plain_2   =  plot(sol_plain.x, sol_plain.y(2,:), 'k.-', 'DisplayName', 'plain ode15s');
+ylabel('i_C (A) = i_L (A) (via constraint)');
+Switch_plot_2  = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+xlabel('Time (s)');
+legend();
+hold off;
+```
+
+![](plots_rlc/rlc_plot_1.png)
+![](plots_rlc/rlc_plot_close.png)
+
+To further investigate this example, take a look at the files `rlc_main.m` and `rlcRHs.m`.
+
+## Sources
+
+- Physics Foundation: Alexander, C. K., & Sadiku, M. N. O. Fundamentals of Electric Circuits.
+- DAE Formulation: Kuzmenko, D. (2018). Switched nonlinear DAEs in electrical circuit theory. Bachelor Thesis, supervised by Prof. Dr. Stephan Trenn.

toolbox/examples/rlcExample/rlcRHS.m

@@ -0,0 +1,14 @@
+function dx = rlcRHS(~,x,p)
+    dx  = zeros(3,1);
+    L   = p(1); R1  = p(2); R2  = p(3);
+    C   = p(4); Vs  = p(5); Vth = p(6);
+    iL  = x(1); VC  = x(2); iC  = x(3);
+
+    if VC > Vth
+        dx(1) = (Vs - R1*iL - VC)/L;
+    else
+        dx(1) = (Vs - R2*iL - VC)/L;
+    end 
+    dx(2) = iC/C;
+    dx(3) = iL - iC; % algebraic constraint
+end

toolbox/examples/rlcExample/rlc_main.m

@@ -0,0 +1,48 @@
+% RLC circuit with fuse
+% for more details, see rlcExample_README.md
+
+%% Setup and integration
+integrator = @ode15s;
+M          = diag([1,1,0]);
+x0         = [0; 0; 0];
+tspan      = [0 5]; 
+
+L   = 1.0;
+R1  = 5.0;
+R2  = 0.1;
+C   = 1.0;
+Vs  = 1.0;
+Vth = 0.5;
+p   = [L; R1; R2; C; Vs; Vth];
+
+opts_ifdiff = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5);
+opts_plain  = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5);
+
+datahandle = prepareDatahandleForIntegration('rlcRHS', 'integrator', integrator, 'options', opts_ifdiff);
+
+sol_ifdiff = solveODE(datahandle, tspan, x0, p);
+sol_plain  = integrator(@(t, x) rlcRHS(t, x, p), tspan, x0, opts_plain);
+
+%% Plots
+clf;
+fig1 = figure(01);
+
+subplot(2,1,1);
+hold on;
+Plot_ifdiff_1   = plot(sol_ifdiff.x, sol_ifdiff.y(1,:), 'ro--', 'DisplayName', 'IFDIFF'); 
+Plot_plain_1    = plot(sol_plain.x, sol_plain.y(1,:), 'k.-', 'DisplayName', 'plain ode15s');
+%Switch_plot_1   = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+ylabel('i_L (A)');
+xlabel('Time (s)');
+legend();
+hold off;
+
+subplot(2,1,2);
+hold on
+Plot_ifdiff_2  = plot(sol_ifdiff.x, sol_ifdiff.y(2,:), 'ro--', 'DisplayName', 'IFDIFF' );
+Plot_plain_2   =  plot(sol_plain.x, sol_plain.y(2,:), 'k.-', 'DisplayName', 'plain ode15s');
+ylabel('i_C (A) = i_L (A) (via constraint)');
+Switch_plot_2  = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch');
+xlabel('Time (s)');
+legend();
+hold off;

toolbox/internal/solving/extendODEuntilSwitch.m

@@ -35,6 +35,9 @@ function extendODEuntilSwitch(datahandle)
 t2FromRootFinding = data.SWP_detection.t2;
 baseOffset = 16*eps(data.SWP_detection.t2);
 iter = 0;
+
+%disp(config.switchingPointErrorThreshold)
+
 while isempty(switchingIndices)
     data = datahandle.getData();
 
@@ -51,6 +54,7 @@ function extendODEuntilSwitch(datahandle)
     % We can not start integrating from the old t2 because this would result in integration over a tiny
     % interval whose result would vanish due to limited floating point accuracy.
     data.SWP_detection.t2 = data.SWP_detection.t2 + baseOffset * 10^min(iter, config.switchingPointMaxPower);
+
     datahandle.setData(data);
 
     % Check if there is a new signature.