diff --git a/README.md b/README.md index d0ee00e..1c7f27d 100644 --- a/README.md +++ b/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); ``` diff --git a/toolbox/examples/daeExample/daeExampleRHS.m b/toolbox/examples/daeExample/daeExampleRHS.m new file mode 100644 index 0000000..a13027e --- /dev/null +++ b/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 diff --git a/toolbox/examples/daeExample/daeExample_README.md b/toolbox/examples/daeExample/daeExample_README.md new file mode 100644 index 0000000..821c500 --- /dev/null +++ b/toolbox/examples/daeExample/daeExample_README.md @@ -0,0 +1,112 @@ + + + +# 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. diff --git a/toolbox/examples/daeExample/daeExample_main.m b/toolbox/examples/daeExample/daeExample_main.m new file mode 100644 index 0000000..98dacc5 --- /dev/null +++ b/toolbox/examples/daeExample/daeExample_main.m @@ -0,0 +1,25 @@ +% DAE Example Main + +%% 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; diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1.png b/toolbox/examples/daeExample/plots_daeExample/plot1.png new file mode 100644 index 0000000..e88b597 Binary files /dev/null and b/toolbox/examples/daeExample/plots_daeExample/plot1.png differ diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1_2.png b/toolbox/examples/daeExample/plots_daeExample/plot1_2.png new file mode 100644 index 0000000..5b65729 Binary files /dev/null and b/toolbox/examples/daeExample/plots_daeExample/plot1_2.png differ diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1_close.png b/toolbox/examples/daeExample/plots_daeExample/plot1_close.png new file mode 100644 index 0000000..fa5aa7b Binary files /dev/null and b/toolbox/examples/daeExample/plots_daeExample/plot1_close.png differ diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1_close2.png b/toolbox/examples/daeExample/plots_daeExample/plot1_close2.png new file mode 100644 index 0000000..0f9ce45 Binary files /dev/null and b/toolbox/examples/daeExample/plots_daeExample/plot1_close2.png differ diff --git a/toolbox/examples/daeExample/plots_daeExample/switching_point_ifdiff_plain.fig b/toolbox/examples/daeExample/plots_daeExample/switching_point_ifdiff_plain.fig new file mode 100644 index 0000000..9ad10f3 Binary files /dev/null and b/toolbox/examples/daeExample/plots_daeExample/switching_point_ifdiff_plain.fig differ diff --git a/toolbox/examples/rlcExample/plots_rlc/plot_rlc.png b/toolbox/examples/rlcExample/plots_rlc/plot_rlc.png new file mode 100644 index 0000000..5486a71 Binary files /dev/null and b/toolbox/examples/rlcExample/plots_rlc/plot_rlc.png differ diff --git a/toolbox/examples/rlcExample/plots_rlc/rlc_plot_1.png b/toolbox/examples/rlcExample/plots_rlc/rlc_plot_1.png new file mode 100644 index 0000000..d6b6f8d Binary files /dev/null and b/toolbox/examples/rlcExample/plots_rlc/rlc_plot_1.png differ diff --git a/toolbox/examples/rlcExample/plots_rlc/rlc_plot_close.png b/toolbox/examples/rlcExample/plots_rlc/rlc_plot_close.png new file mode 100644 index 0000000..bba7f40 Binary files /dev/null and b/toolbox/examples/rlcExample/plots_rlc/rlc_plot_close.png differ diff --git a/toolbox/examples/rlcExample/rlcExample_README.md b/toolbox/examples/rlcExample/rlcExample_README.md new file mode 100644 index 0000000..74df167 --- /dev/null +++ b/toolbox/examples/rlcExample/rlcExample_README.md @@ -0,0 +1,119 @@ + + +# 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. \ No newline at end of file diff --git a/toolbox/examples/rlcExample/rlcRHS.m b/toolbox/examples/rlcExample/rlcRHS.m new file mode 100644 index 0000000..2aaead8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/toolbox/examples/rlcExample/rlc_main.m b/toolbox/examples/rlcExample/rlc_main.m new file mode 100644 index 0000000..0791195 --- /dev/null +++ b/toolbox/examples/rlcExample/rlc_main.m @@ -0,0 +1,47 @@ +% RLC circuit with fuse + +%% 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; diff --git a/toolbox/internal/solving/extendODEuntilSwitch.m b/toolbox/internal/solving/extendODEuntilSwitch.m index 28d96be..4fc986d 100644 --- a/toolbox/internal/solving/extendODEuntilSwitch.m +++ b/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. diff --git a/toolbox/internal/solving/extendODEuntilSwitch_t1_to_t2.m b/toolbox/internal/solving/extendODEuntilSwitch_t1_to_t2.m index f6fe78b..81891a2 100644 --- a/toolbox/internal/solving/extendODEuntilSwitch_t1_to_t2.m +++ b/toolbox/internal/solving/extendODEuntilSwitch_t1_to_t2.m @@ -6,23 +6,46 @@ function extendODEuntilSwitch_t1_to_t2(datahandle) config = makeConfig(); data = datahandle.getData(); -t = data.SWP_detection.solution_until_t1.x(end); -x = deval(data.SWP_detection.solution_until_t1, data.SWP_detection.solution_until_t1.x(end)); +% last integrator step before switch +ti = data.SWP_detection.solution_until_t1.x(end); +% solution at the last time point before switch +x = deval(data.SWP_detection.solution_until_t1, ti); solution = data.SWP_detection.solution_until_t1; + end_point = data.SWP_detection.t2; -options = data.integratorSettings.options; -ctrlif_setForcedBranchingSignature(datahandle, t, x); +% new step size +delta_t = end_point - ti; + +options = data.integratorSettings.options; +solver = solution.solver; + +% Common setup for all solvers +ctrlif_setForcedBranchingSignature(datahandle, ti, x); data = datahandle.getData(); data.caseCtrlif = config.caseCtrlif.extendODEuntilSwitch; datahandle.setData(data); z = odextend(solution, [], end_point, [], options); - data = datahandle.getData(); data.SWP_detection.solution_until_t2 = z; +% last point strategy for one step solvers +compatible_solvers = {'ode23', 'ode45', 'ode78', 'ode89', 'ode113', 'ode23t', 'ode23tb'}; +if config.last_point_strategy && ismember(solver, compatible_solvers) + + data.integratorSettings.options.InitialStep = delta_t; + data.integratorSettings.options.AbsTol = 1; + data.integratorSettings.options.RelTol = 1; + + datahandle.setData(data); + + z = odextend(solution, [], end_point, [], options); + data = datahandle.getData(); + data.SWP_detection.solution_until_t2 = z; +end + datahandle.setData(data); end diff --git a/toolbox/makeConfig.m b/toolbox/makeConfig.m index d235daf..d4019c9 100644 --- a/toolbox/makeConfig.m +++ b/toolbox/makeConfig.m @@ -139,5 +139,8 @@ config.swfreq_haltOnWarning = false; +% last point strategy for one-step solvers +config.last_point_strategy = true; + config_out = config; end