diff --git a/toolbox/examples/predatorpreyFilippov/pprhs.m b/toolbox/examples/predatorpreyFilippov/pprhs.m deleted file mode 100644 index e1b1a8ca..00000000 --- a/toolbox/examples/predatorpreyFilippov/pprhs.m +++ /dev/null @@ -1,44 +0,0 @@ -function dx = pprhs(t,x,p) -% Predator Prey Model -% -% Source: -% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: -% Sliding Shilnikov connection in Filippov-type predator-prey model -% Nonlinear Dyn (2020) 100:2973-2987 -% - -% preallocate output -dx = [0 ; 0 ; 0]; - -% x = (p1, p2, P) % prey, prey2, predator -p1 = x(1); -p2 = x(2); -P = x(3); - -% p = (r1, r2, beta1, beta2, q1, q2, m, e, aq) % parameters -r1 = p(1); -r2 = p(2); -beta1 = p(3); -beta2 = p(4); -q1 = p(5); -q2 = p(6); -m = p(7); -e = p(8); -aq = p(9); - -% switching manifold -H = beta1*p1 - aq*beta2*p2; - -% switched rhs -if H >= 0 - dx(1) = (r1 - beta1*P) .* p1; - dx(2) = r2*p2; - dx(3) = (e*q1*beta1*p1 - m) .* P; -else - dx(1) = r1*p1; - dx(2) = (r2 - beta2*P) .* p2; - dx(3) = (e*q2*beta2*p2 - m) .* P; -end - - -end \ No newline at end of file diff --git a/toolbox/examples/predatorpreyFilippov/pprhs5.m b/toolbox/examples/predatorpreyFilippov/pprhs5.m deleted file mode 100644 index b44d418b..00000000 --- a/toolbox/examples/predatorpreyFilippov/pprhs5.m +++ /dev/null @@ -1,49 +0,0 @@ -function dx = pprhs5(t,state,p) -% Predator Prey Model -% -% Source: -% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: -% Sliding Shilnikov connection in Filippov-type predator-prey model -% Nonlinear Dyn (2020) 100:2973-2987 -% -% FORMULA 5 -% X(x,y,z) = [ (r1-z)*x , r2*y , (e*q1*x-m)*z ]; -% Y(x,y,z) = [ r1*x , (r2-beta2/beta1*z)*y , (e*q2/aq*y-m)*z ]; -% - -% preallocate output -dx = [0 , 0 , 0]; - -% x = (p1, p2, P) % prey, prey2, predator -x = state(1); -y = state(2); -z = state(3); - -% p = (r1, r2, beta1, beta2, q1, q2, m, e, aq) % parameters -r1 = p(1); -r2 = p(2); -beta1 = p(3); -beta2 = p(4); -q1 = p(5); -q2 = p(6); -m = p(7); -e = p(8); -aq = p(9); - -% switching manifold -h = x - y; - -% switched rhs -if h > 0 - % X(x,y,z) = [ (r1-z)*x , r2*y , (e*q1*x-m)*z ]; - dx(1) = (r1-z)*x ; - dx(2) = r2*y ; - dx(3) = (e*q1*x-m)*z ; -else - % Y(x,y,z) = [ r1*x , (r2-beta2/beta1*z)*y , (e*q2/aq*y-m)*z ]; - dx(1) = r1*x ; - dx(2) = (r2-beta2/beta1*z)*y ; - dx(3) = (e*q2/aq*y-m)*z ; -end - -end \ No newline at end of file diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_formula5_rhs.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_formula5_rhs.m new file mode 100644 index 00000000..f4aca2ca --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_formula5_rhs.m @@ -0,0 +1,45 @@ +function dx = predatorPrey3D_formula5_rhs(~, state, p) +% Predator Prey Model +% +% Source: +% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: +% Sliding Shilnikov connection in Filippov-type predator-prey model +% Nonlinear Dyn (2020) 100:2973-2987 +% +% FORMULA 5 +% X(x,y,z) = [ (r1-z)*x , r2*y , (e*q1*x-m)*z ]; +% Y(x,y,z) = [ r1*x , (r2-beta2/beta1*z)*y , (e*q2/aq*y-m)*z ]; + +% Preallocate +dx = zeros(3, 1); + +% Unpack state components and parameters (notation from paper) +x = state(1); % prey1 +y = state(2); % prey2 +z = state(3); % predator + +r1 = p(1); +r2 = p(2); +beta1 = p(3); +beta2 = p(4); +q1 = p(5); +q2 = p(6); +m = p(7); +e = p(8); +aq = p(9); + +% Switching manifold +H = x - y; + +if H > 0 + % X(x,y,z) = [ (r1-z)*x , r2*y , (e*q1*x-m)*z ]; + dx(1) = (r1-z) * x; + dx(2) = r2 * y; + dx(3) = (e*q1*x - m) * z; +else + % Y(x,y,z) = [ r1*x , (r2-beta2/beta1*z)*y , (e*q2/aq*y-m)*z ]; + dx(1) = r1 * x; + dx(2) = (r2 - beta2/beta1*z) * y; + dx(3) = (e*q2/aq*y - m) * z; +end +end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m new file mode 100644 index 00000000..68f3fefd --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m @@ -0,0 +1,39 @@ +function dx = predatorPrey3D_rhs(~, x, p) +% 2-Prey-1-Predatory Model +% +% Source: +% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: +% Sliding Shilnikov connection in Filippov-type predator-prey model +% Nonlinear Dyn (2020) 100:2973-2987 + +% Preallocate +dx = zeros(3, 1); + +% Unpack state components and parameters +p1 = x(1); % prey1 +p2 = x(2); % prey2 +P = x(3); % Predator + +r1 = p(1); +r2 = p(2); +beta1 = p(3); +beta2 = p(4); +q1 = p(5); +q2 = p(6); +m = p(7); +e = p(8); +aq = p(9); + +% Switching manifold +H = beta1*p1 - aq*beta2*p2; + +if H >= 0 + dx(1) = (r1 - beta1*P) .* p1; + dx(2) = r2*p2; + dx(3) = (e*q1*beta1*p1 - m) .* P; +else + dx(1) = r1*p1; + dx(2) = (r2 - beta2*P) .* p2; + dx(3) = (e*q2*beta2*p2 - m) .* P; +end +end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m new file mode 100644 index 00000000..20d7a2fd --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m @@ -0,0 +1,217 @@ +%% Setup +% Time span +tspan = [0, 100]; + +% Integrator +intIfdiff = @ode45; +intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12); +eulerStep = 1e-7; + +% Parameter values required for Shilnikov behavior in the paper (see RHS file). +m = 0.790; +r1 = 0.836; +e = 0.948; +q1 = 0.772; +aq = 0.660; +beta2 = 0.896; +% Some values changed for prettier plot (remains a Filippov system). +q2 = 1.500; % 1.084 (paper) +r2 = 0.300; % 0.126 (paper) +beta1 = 7.810; % 7.800 (Example plot in paper, but can be chosen arbitrarily, retaining Shilnikov behavior) +p = [r1, r2, beta1, beta2, q1, q2, m, e, aq]; + +% Initial values (used for plot in paper, but may be changed) +a = 0.286975; +x0 = [a; a; r1-r2]; + +% RHS function dx = f(t, x, p), must be implemented in separate file for IFDIFF. +rhs = @predatorPrey3D_rhs; + +% Plotting +plot_n = 10000; +plot_t = linspace(tspan(1), tspan(end), plot_n); + +nameIfdiff = sprintf('Ifdiff/%s', func2str(intIfdiff)); +nameEuler = sprintf('Euler/%.0e', eulerStep); +lwIfdiff = 3; +lwEuler = 3; +colorIfdiff = [0, 0.447, 0.741]; +colorEuler = [0.85, 0.325, 0.098]; +lsIfdiff = '-'; +lsEuler = '--'; + +% Sensitivity +wrt_y = 3; +eulerDisturbH = 1e-6; + + +%% Solve with IFDIFF +datahandle = prepareDatahandleForIntegration(rhs, 'solver', intIfdiff, 'options', intOptions); +% Enable storing of convexification data +configNew = makeConfig(); +configNew.storeSlidingInfo = true; +configOld = makeConfig(configNew); +try + solIfdiff = solveODE(datahandle, tspan, x0, p); +catch ME + makeConfig(configOld); + rethrow(ME); +end +makeConfig(configOld); + +%% Solve with Euler +[owndir, ~] = fileparts(mfilename('fullpath')); +if contains(owndir, 'Editor_') + [owndir, ~] = fileparts(matlab.desktop.editor.getActiveFilename); +end + +euler_prefix = sprintf('sol_euler_%.0e', eulerStep); +euler_vname = 'sol_euler'; +euler_fname = fullfile(owndir, [euler_prefix, '.mat']); + +solEuler = loadEulerOrCompute(euler_fname, euler_vname, rhs, tspan, x0, p, eulerStep); + + +%% Plot solution +plot_x_ifdiff = deval(solIfdiff, plot_t); +axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff); + +plot_x_euler = interp1(solEuler.x, solEuler.y', plot_t)'; +plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler); + + +%% Compute Sensitivity with IFDIFF +FDstep = generateFDstep(length(x0), length(p)); +sensFunVDE = generateSensitivityFunction(datahandle, solIfdiff, FDstep, ... + 'method', 'VDE', ... + 'calcGy', true, ... + 'calcGp', false); +sensVDE = sensFunVDE(plot_t); +G = arrayfun(@(s) s.Gy(:, wrt_y), sensVDE, 'UniformOutput', false); +plot_sens_ifdiff = [G{:}]; + + +%% Compute Sensitivity with Euler +euler_prefix_disturb = sprintf('%s_disturb_wrt_y%d', euler_prefix, wrt_y); +euler_fname_disturb = fullfile(owndir, [euler_prefix_disturb, '.mat']); + +x0_disturb = x0; +x0_disturb(wrt_y) = x0_disturb(wrt_y) + eulerDisturbH; +solEulerDisturb = loadEulerOrCompute(euler_fname_disturb, euler_vname, rhs, tspan, x0_disturb, p, eulerStep); +sensEuler = (solEulerDisturb.y - solEuler.y) / eulerDisturbH; +plot_sens_euler = interp1(solEuler.x, sensEuler', plot_t)'; + + +%% Plot Sensitivity +for idx_y=1:3 +ax = plotSens([], plot_t, plot_sens_ifdiff(idx_y, :), 'sensIFDIFF', lwIfdiff, colorIfdiff, lsIfdiff); +plotSens(ax, plot_t, plot_sens_euler(idx_y, :), 'sensEuler', lwEuler, colorEuler, lsEuler); +end + + +%% Plot alpha and switching function +data = datahandle.getData(); +t_alpha = data.sliding.convexification.t; +alpha = data.sliding.convexification.alpha; +ax = plotSwitchingFuncOrAlpha([], t_alpha, alpha, 'Alpha', 2, colorIfdiff, 'left', 'Alpha'); + +sw = solIfdiff.switchingFunction{1}; +sw_x = arrayfun(@(t) sw([], t, deval(solIfdiff, t), p), t_alpha); +plotSwitchingFuncOrAlpha(ax, t_alpha, sw_x, 'Switching Function', 2, colorEuler, 'right', 'Switching Function'); + +plotSwitches(ax, solIfdiff.switches) + + +%% Helpers +function sol_euler = loadEulerOrCompute(fname, vname, rhs, tspan, x0, p, step) +if isfile(fname) + tmp = load(fname, vname); + sol_euler = tmp.(vname); + fprintf('Loading %s from file %s\n', vname, fname); + return +end +fprintf('Integrating with explicit Euler (might take a while) ...\n') +sol_euler = explEuler(rhs, tspan, x0, p, step); +fprintf('Saving result to %s for later reuse.\n', fname); +save(fname, vname); +end + +function plotSwitches(ax, switches) +xline(ax, switches, '--', 'LineWidth', 1.5, 'HandleVisibility', 'off'); +end + +function ax = plotSwitchingFuncOrAlpha(ax, t, x, name, lw, color, yyax, ylab) +if isempty(ax) + f = figure; + ax = axes(f); +end +yyaxis(ax, yyax); +hold(ax, 'on'); +plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color); +hold(ax, 'off'); +grid(ax, 'on'); +xlabel(ax, 't'); +ylabel(ax, ylab); +legend(ax, 'location', 'northeast'); + +yl = ylim(ax); +yl(1) = 0; +ylim(ax, yl); +end + +function ax = plotSens(ax, t, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +hold(ax, 'on'); +plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls); +hold(ax, 'off'); +grid(ax, 'on'); +xlabel(ax, 'Time') +ylabel(ax, 'Sensitivity') +legend(ax, 'location', 'northeast'); +end + +function ax = plotSol3d(ax, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +hold(ax, 'on'); +plot3(ax, x(3, :), x(2, :), x(1, :), 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls); +hold(ax, 'off'); + +view(ax, [97 51]); +box(ax, 'on'); +grid(ax, 'on'); +xlabel(ax, 'Predator'); +ylabel(ax, 'Prey 2'); +zlabel(ax, 'Prey 1'); +legend(ax, 'location', 'northeast'); +end + +function sol = explEuler(rhs, tspan, x0, p, stepsize) +xdim = length(x0); +stepcount = (tspan(end)-tspan(1)) / stepsize; +sfac = 0.001; % store factor +n_out = ceil(stepcount*sfac) + 1; + +Xi = reshape(x0, [], 1); +X = zeros(xdim, n_out); +X(:,1) = Xi; + +k = 2; nextout = ceil(1 / sfac); +for i=2:stepcount + Xi = Xi + stepsize * rhs(i*stepsize, Xi, p); + if (i == nextout) + X(:,k) = Xi; k = k + 1; + nextout = nextout + ceil(1 / sfac); + end + if ~mod(i, stepcount / 100), fprintf('.'); end +end +fprintf('\n'); +T = linspace(tspan(1), tspan(end), n_out); +sol.x = T; +sol.y = X; +end diff --git a/toolbox/examples/predatorpreyFilippov/testpprhs.m b/toolbox/examples/predatorpreyFilippov/testpprhs.m deleted file mode 100644 index 1dc27952..00000000 --- a/toolbox/examples/predatorpreyFilippov/testpprhs.m +++ /dev/null @@ -1,162 +0,0 @@ -%% SETUP -% parameter values p = (r1, r2, beta1, beta2, q1, q2, m, e) -% order as in paper (see rhs file) -m = 0.790; -r1 = 0.836; -e = 0.948; -q1 = 0.772; -aq = 0.660; -beta2 = 0.896; - -beta1 = 7.81; q2 = 1.5; r2 = 0.3; - -% initial values, parameters, timespan -a = 0.286975; -x0_1 = [a ; a ; r1-r2]; -p = [r1, r2, beta1, beta2, q1, q2, m, e, aq]; -tspan = [0 100]; - -% configure plotting -X_plot = linspace(tspan(1), tspan(end), 1000); -fignum = 1000; -figure(fignum); clf; hold('on'); -plotit = @plotter; - -% solver selection and configuration -intEuler = @explEuler; -eulerStep = 1e-7; -intIfdiff = @ode45; -intOptions = odeset('reltol', 1e-5, 'abstol', 1e-5, 'MaxStep', 0.5); - -% select what to do -doIfdiff = true; -doEuler = false; % Euler takes very long! false: loads from file, if exists -doErrorplot = true; % compare ifdiff and plain euler - -% name generators -nameIfdiff = @(f) sprintf('ifdiff/%s', func2str(f)); -namePlain = @(f) sprintf('plain %s' , func2str(f)); - - -%% IFDIFF -if doIfdiff - fprintf('Integrating with IFDIFF/%s ...\n', func2str(intIfdiff)) - figure(fignum); - datahandle = prepareDatahandleForIntegration('pprhs', 'solver', func2str(intIfdiff), 'options', intOptions); - th = tic(); - sol_ifdiff = solveODE(datahandle, tspan, x0_1, p); - time_ifdiff = toc(th); fprintf('IFDIFF took %g s\n', time_ifdiff); - X_ifdiff = X_plot; - Y_ifdiff = deval(sol_ifdiff, X_ifdiff); - linewidth = 3.0; - hIFDIFF = plotit(fignum, Y_ifdiff, 'g', nameIfdiff(intIfdiff), linewidth); -end - - -%% EULER Integration -[owndir, ~] = fileparts(mfilename('fullpath')); -euler_fname = fullfile(owndir, sprintf('sol_euler_%.0e.mat', eulerStep)); -if doEuler - fprintf('Integrating with integrator %s (might take a while) ...\n', func2str(intEuler)) - figure(fignum); - th = tic(); - sol_euler = intEuler(@(t,x) pprhs(t,x,p), tspan, x0_1, eulerStep); - time_euler = toc(th); fprintf('Euler took %g s\n', time_euler); - fprintf('Saving result to %s for later reuse.\n', euler_fname); - save(euler_fname, "sol_euler"); -elseif isfile(euler_fname) - fprintf('Loading sol_euler from file %s\n', euler_fname); - tmp = load(euler_fname, 'sol_euler'); - sol_euler = tmp.sol_euler; - doEuler = true; -end -if doEuler - X_euler = X_plot; - Y_euler = transpose(interp1(sol_euler.x, transpose(sol_euler.y), X_euler)); - linewidth = 2.0; - hEuler = plotit(fignum, Y_euler, 'c', namePlain(intEuler), linewidth); -end - -%% ANALYSIS -if doErrorplot - if (doIfdiff && doEuler) - fignum = fignum + 1; errorPlot(fignum, X_plot, Y_ifdiff, Y_euler, nameIfdiff(intIfdiff), namePlain(intEuler)) - else - fprintf('Not enough data to compare. Either IFDIFF or Euler solution missing.\n') - end -end - -% FINITO -return - - -%% HELPERS -function errorPlot(fignum, x1, y1, y2, intname1, intname2) - figure(fignum); clf(fignum); - ydiff = calcDiff(y1, y2); - semilogy(x1, ydiff, 'LineWidth', 1.0); - xlabel('t'); - ylabel('||y||_2'); - title(sprintf('difference %s and %s', intname1, intname2)); - drawnow -end - - -function h = plotter(fignum, y, color, name, lw) - figure(fignum); hold on; - h = plot3(y(3,:), y(2,:), y(1,:), 'Color', color, 'LineWidth', lw, 'DisplayName', name); - view([97 51]); - grid on; - box on; - xlabel('Predator'); - ylabel('Prey 2'); - zlabel('Prey 1'); - legend('location', 'northeast'); - drawnow - pause(1.0); - set(fignum, 'Position', [200 250 750 375]); -end - - -function diffnorm = calcDiff(yA, yB) - % Go through x-coordinates of yA, find nearest x-coordinate in yB and compare. - % Algorithm: 1) Go through x of yA in Index i - % 2) Find nearest x in yB in a window centered around yB(:,i) - len = length(yA); - diffnorm = zeros(len, 1); - for i = 1:len - y = yA(:, i); - % find nearest x in yB in a time-window around the timepoint of y - window = 250; - j0 = max(1, i-window); - jf = max(len, j0+window); - jidx = j0:jf; - tmpdiff = yB(:, jidx) - repmat(y, [1, length(jidx)]); - tmpdiff = vecnorm(tmpdiff,2,1); - diffnorm(i) = min(tmpdiff); - diffnorm(i) = diffnorm(i) / norm(y); - end -end - - -function sol = explEuler(rhs, tspan, x0, stepsize) - xdim = length(x0); % get dimension - stepcount = (tspan(end)-tspan(1))/stepsize; - sfac = 0.001; % store factor - X = zeros(xdim, ceil(stepcount*sfac)+1); - Xi = reshape(x0, [], 1); - X(:,1) = Xi; - k = 2; nextout = ceil(1 / sfac); - for i=2:stepcount - Xi = Xi + stepsize * rhs(i*stepsize, Xi); - if (i == nextout) - X(:,k) = Xi; k = k + 1; - nextout = nextout + ceil(1 / sfac); - end - if ~mod(floor(100*i/stepcount), 10), fprintf('.'); end - end - fprintf('\n') - T = linspace(tspan(1), tspan(end), ceil(stepcount*sfac)+1); - sol.x = T; - sol.y = X; -end diff --git a/toolbox/internal/functionGeneration/BranchingSignature.m b/toolbox/internal/functionGeneration/BranchingSignature.m index 2eb0e2ad..fa63a065 100644 --- a/toolbox/internal/functionGeneration/BranchingSignature.m +++ b/toolbox/internal/functionGeneration/BranchingSignature.m @@ -121,6 +121,7 @@ end end + %% Getters and Setters % Getters retrieve the cached value for str and hash. function str = get.str(this) str = this.strCached; @@ -150,5 +151,10 @@ this.functionIndex = val; this = this.updateCache(); end + + %% Comparison operators + % For equality checks use isequal (checks for matching classes and properties, also works for arrays). + % NOTE: Checked properties seem to include private and hidden properties. + % This should not be an issue, since here the public properties matching ensures the private ones matching also. end end diff --git a/toolbox/internal/sensitivities/getGy_Gp_update.m b/toolbox/internal/sensitivities/getGy_Gp_update.m index 838ed61b..c053ff5b 100644 --- a/toolbox/internal/sensitivities/getGy_Gp_update.m +++ b/toolbox/internal/sensitivities/getGy_Gp_update.m @@ -1,91 +1,93 @@ function Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) - % Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) - % - % Calculates the update-matrices for the sensitivitiy calculation with END_piecewise or VDE. - % - % INPUT: datahandle - datahandle you get from the integration process with solveODE - % amountG - amount of intermediate G-matrices that already have been calculated and saved from previous function calls. - % (If amountG = 1, it means that the intermediates were only initialised and all the matrices have still to be calculated.) - % modelNum - number of the model until which the update-matrices need to be calculated - % Gp_flag - flag that is true if the sensitivities with respect to the parameters should be calculated - % options - struct that has the informations for the stepsizes for END, the integrator and the method - % - % OUTPUT: Updates - struct that has he update-matrices needed for the sensitivities w.r.t. y0 and optionally w.r.t. p +% Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) +% +% Calculates the update-matrices for the sensitivitiy calculation with END_piecewise or VDE. +% +% INPUT: datahandle - datahandle you get from the integration process with solveODE +% amountG - amount of intermediate G-matrices that already have been calculated and saved from previous function calls. +% (If amountG = 1, it means that the intermediates were only initialised and all the matrices have still to be calculated.) +% modelNum - number of the model until which the update-matrices need to be calculated +% Gp_flag - flag that is true if the sensitivities with respect to the parameters should be calculated +% options - struct that has the informations for the stepsizes for END, the integrator and the method +% +% OUTPUT: Updates - struct that has he update-matrices needed for the sensitivities w.r.t. y0 and optionally w.r.t. p - % Initialize - data = datahandle.getData(); - parameters = data.SWP_detection.parameters; - switches = data.computeSensitivity.switches_extended; - switches_left = data.computeSensitivity.switches_extended_left; - y_to_switches = data.computeSensitivity.y_to_switches; - y_to_switches_left = data.computeSensitivity.y_to_switches_left; - switching_functions = data.SWP_detection.switchingFunction; - jump_functions = data.SWP_detection.jumpFunction; - functionRHS = data.integratorSettings.preprocessed_rhs; +% Initialize +data = datahandle.getData(); +parameters = data.SWP_detection.parameters; +switches = data.computeSensitivity.switches_extended; +switches_left = data.computeSensitivity.switches_extended_left; +y_to_switches = data.computeSensitivity.y_to_switches; +y_to_switches_left = data.computeSensitivity.y_to_switches_left; +switching_functions = data.SWP_detection.switchingFunction; +jump_functions = data.SWP_detection.jumpFunction; - FDstep = options.FDstep; +FDstep = options.FDstep; - dim_y = data.computeSensitivity.dim_y; - dim_p = data.computeSensitivity.dim_p; - unit = eye(dim_y); +dim_y = data.computeSensitivity.dim_y; +dim_p = data.computeSensitivity.dim_p; +unit = eye(dim_y); - numModels = endModel - startModel; - Updates.Uy_new = cell(1, numModels); +numModels = endModel - startModel; +Updates.Uy_new = cell(1, numModels); - if Gp_flag - Updates.Up_new = cell(1, numModels); - h_p = fdStep_getH_p(FDstep, parameters); - end - - % Fix the model and set the model number for which model you want to evaluate the RHS - data.computeSensitivity.modelStage = startModel; - datahandle.setData(data); +if Gp_flag + Updates.Up_new = cell(1, numModels); + h_p = fdStep_getH_p(FDstep, parameters); +end - % Calculate the remaining update matrices for the update formula until modelNum. - for i = startModel : (endModel-1) - switchingFunction = switching_functions{i}; - jumpFunction = jump_functions{i}; - ts = switches(i+1); - tsMinus = switches_left(i+1); - yMinus = y_to_switches_left(:, i+1); - yPlus = y_to_switches(:, i+1); +% Fix the model and set the model number for which model you want to evaluate the RHS +data.computeSensitivity.modelStage = startModel; +datahandle.setData(data); - h_y = fdStep_getH_y(FDstep, yMinus); - h_t = fdStep_getH_t(FDstep, tsMinus); +rhs = getRhsFromModelNum(datahandle, startModel); +% Calculate the remaining update matrices for the update formula until modelNum. +for i = startModel : (endModel-1) + switchingFunction = switching_functions{i}; + jumpFunction = jump_functions{i}; + ts = switches(i+1); + tsMinus = switches_left(i+1); + yMinus = y_to_switches_left(:, i+1); + yPlus = y_to_switches(:, i+1); - % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) - del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); - del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); - diff_sigmat = del_sigmat + del_sigmay * functionRHS(datahandle, tsMinus, yMinus, parameters); - if isempty(jumpFunction) - del_jumpy = zeros(dim_y); - del_jumpt = zeros(dim_y, 1); - else - del_jumpy = del_f_del_y(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_y); - del_jumpt = del_f_del_t(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_t); - end + h_y = fdStep_getH_y(FDstep, yMinus); + h_t = fdStep_getH_t(FDstep, tsMinus); - % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number - % and evalate the RHS with the model fixed on the left of the switch. - fminus = functionRHS(datahandle, tsMinus, yMinus, parameters); + % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) + del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); + del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); + diff_sigmat = del_sigmat + del_sigmay * rhs(datahandle, tsMinus, yMinus, parameters); + if isempty(jumpFunction) + del_jumpy = zeros(dim_y); + del_jumpt = zeros(dim_y, 1); + else + del_jumpy = del_f_del_y(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_y); + del_jumpt = del_f_del_t(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_t); + end - data = datahandle.getData(); - data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; - datahandle.setData(data); + % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number + % and evalate the RHS with the model fixed on the left of the switch. + fminus = rhs(datahandle, tsMinus, yMinus, parameters); - fplus = functionRHS(datahandle, ts, yPlus, parameters); + % If the RHS has ctrlifs, this ensures that the correct branching is applied. + data = datahandle.getData(); + data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; + datahandle.setData(data); + % Get RHS for next submodel (will correspond to fminus in the next iteration) + rhs = getRhsFromModelNum(datahandle, i+1); + fplus = rhs(datahandle, ts, yPlus, parameters); - % Calculate the updates according to the update formula - Updates.Uy_new{i} = unit + del_jumpy + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmay / diff_sigmat; + % Calculate the updates according to the update formula + Updates.Uy_new{i} = unit + del_jumpy + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmay / diff_sigmat; - if Gp_flag - del_sigmap = del_f_del_p(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_p); - if isempty(jumpFunction) - del_jumpp = zeros(dim_y, dim_p); - else - del_jumpp = del_f_del_p(datahandle, jumpFunction, tsMinus, yMinus, parameters, h_p); - end - Updates.Up_new{i} = del_jumpp + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmap / diff_sigmat; + if Gp_flag + del_sigmap = del_f_del_p(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_p); + if isempty(jumpFunction) + del_jumpp = zeros(dim_y, dim_p); + else + del_jumpp = del_f_del_p(datahandle, jumpFunction, tsMinus, yMinus, parameters, h_p); end + Updates.Up_new{i} = del_jumpp + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmap / diff_sigmat; end -end \ No newline at end of file +end +end diff --git a/toolbox/internal/sensitivities/getModelFunction.m b/toolbox/internal/sensitivities/getModelFunction.m deleted file mode 100644 index 8233ccfd..00000000 --- a/toolbox/internal/sensitivities/getModelFunction.m +++ /dev/null @@ -1,35 +0,0 @@ -function functionHandle = getModelFunction(datahandle, modelStage) -%functionHandle = GETMODELFUNCTION(datahandle, modelStage) -% -%Find an existing model function or create a new one based on a model stage. -% -%INPUT: -% datahandle - Contains signature corresponding to the model stage. -% struct -% -% modelStage - Number of the model as encountered during initial computation of solution. -% positive integer -% -%OUTPUT: -% functionHandle - Handle to the requested model function. -% function_handle - -data = datahandle.getData(); -factory = data.codeGen.modelFactory; - -rhsName = data.mtreeplus{2, 1}; -dataSignature = data.SWP_detection.signature; -switchCond = dataSignature.switch_cond{modelStage}; -ctrlifIndex = dataSignature.ctrlif_index{modelStage}; -functionIndex = dataSignature.function_index{modelStage}; -signature = BranchingSignature( ... - rhsName, ... - switchCond, ... - ctrlifIndex, ... - functionIndex); - -[functionHandle, collisionIndex] = factory.findExisting(signature); -if isempty(functionHandle) - functionHandle = factory.createNew(signature, collisionIndex); -end -end diff --git a/toolbox/internal/sensitivities/getRhsFromModelNum.m b/toolbox/internal/sensitivities/getRhsFromModelNum.m new file mode 100644 index 00000000..c1abfa42 --- /dev/null +++ b/toolbox/internal/sensitivities/getRhsFromModelNum.m @@ -0,0 +1,47 @@ +function rhs = getRhsFromModelNum(datahandle, modelNum) +%rhs = GETRHSFROMMODELNUM(datahandle, modelNum) +% +%Return handle to RHS for a particular (potentially Filippov) submodel. +% +%If submodel RHS functions are configured to be exported, this function will attempt +%to find an existing model function or create a new one based on the submodel signature. +% +%INPUT: +% datahandle - Contains stored information about RHS. +% struct +% +% modelNum - Number of the submodel in order of occurrence during initial solution via solveODE. +% positive integer +% +%OUTPUT: +% rhs - Handle to the submodel function corresponding to the signature. +% function_handle + +config = makeConfig(); +data = datahandle.getData(); + +signature = data.SWP_detection.signature{modelNum}; +% Check for Filippov case (multiple signatures) +if ~isscalar(signature) + % TODO: Add option to also export Filippov RHS source code without ctrlifs. + % Assume we don't start in Filippov (guaranteed since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature, ... + switchingFunction, t, y, p); + return +end + +% Get transformed RHS source code without Ctrlifs +if config.removeCtrlifForSensComputation + factory = data.codeGen.modelFactory; + [rhs, collisionIndex] = factory.findExisting(signature); + if isempty(rhs) + rhs = factory.createNew(signature, collisionIndex); + end + return +end + +% Default case, just return the preprocessed RHS +rhs = data.integratorSettings.preprocessed_rhs; diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gp.m b/toolbox/internal/sensitivities/solveDisturbed_Gp.m index 61e2da6e..b3aa3857 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gp.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gp.m @@ -2,32 +2,28 @@ %SOLVEDISTURBED_GP Solve the IVP in the interval tspan with slightly disturbed parameters. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the % parameter disturbed in one component. - config = makeConfig(); - data = datahandle.getData(); - if config.removeCtrlifForSensComputation - functionRHS_original = getModelFunction(datahandle, modelNum); - else - functionRHS_original = data.integratorSettings.preprocessed_rhs; - end - parameters = data.SWP_detection.parameters; - functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); - dim_p = data.computeSensitivity.dim_p; +data = datahandle.getData(); - unit_p = eye(dim_p); +functionRHS_original = getRhsFromModelNum(datahandle, modelNum); - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +parameters = data.SWP_detection.parameters; +functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); +dim_p = data.computeSensitivity.dim_p; - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +unit_p = eye(dim_p); - sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); - sols_disturbed = cell(dim_p, 1); - for j=1:dim_p - parameters_new = parameters + h_p.*unit_p(:,j); - functionRHS_simple_disturb = @(t,y) functionRHS_original(datahandle, t, y, parameters_new); - sols_disturbed{j} = integrator(functionRHS_simple_disturb, tspan, y_start, integrator_options); - end - +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; + +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); + +sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); +sols_disturbed = cell(dim_p, 1); +for j=1:dim_p + parameters_new = parameters + h_p.*unit_p(:,j); + functionRHS_simple_disturb = @(t,y) functionRHS_original(datahandle, t, y, parameters_new); + sols_disturbed{j} = integrator(functionRHS_simple_disturb, tspan, y_start, integrator_options); end +end diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gy.m b/toolbox/internal/sensitivities/solveDisturbed_Gy.m index 48fd2dbd..1a04855b 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gy.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gy.m @@ -2,27 +2,23 @@ %SOLVEDISTURBED_GY Solve the IVP in the interval tspan with slightly disturbed initial values. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the initial % y value disturbed in one component. - config = makeConfig(); - data = datahandle.getData(); - if config.removeCtrlifForSensComputation - functionRHS_original = getModelFunction(datahandle, modelNum); - else - functionRHS_original = data.integratorSettings.preprocessed_rhs; - end - functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); - dim_y = data.computeSensitivity.dim_y; - unit_y = eye(dim_y); +data = datahandle.getData(); - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +functionRHS_original = getRhsFromModelNum(datahandle, modelNum); - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); +dim_y = data.computeSensitivity.dim_y; +unit_y = eye(dim_y); - sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); - sols_disturbed = cell(dim_y, 1); - for j=1:dim_y - sols_disturbed{j} = integrator(functionRHS_simple_END, tspan, y_start + h_y.*unit_y(:,j), integrator_options); - end -end +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; + +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); +sols_disturbed = cell(dim_y, 1); +for j=1:dim_y + sols_disturbed{j} = integrator(functionRHS_simple_END, tspan, y_start + h_y.*unit_y(:,j), integrator_options); +end +end diff --git a/toolbox/internal/sensitivities/solveVDE_Gp.m b/toolbox/internal/sensitivities/solveVDE_Gp.m index 9a656e58..1da1239b 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gp.m +++ b/toolbox/internal/sensitivities/solveVDE_Gp.m @@ -1,25 +1,21 @@ function solVDE = solveVDE_Gp(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GP Solve the VDE for Gp in the interval tspan fixed to model modelNum and return the sol object - config = makeConfig(); - data = datahandle.getData(); +data = datahandle.getData(); - parameters = data.SWP_detection.parameters; - if config.removeCtrlifForSensComputation - rhs = getModelFunction(datahandle, modelNum); - else - rhs = data.integratorSettings.preprocessed_rhs; - end - dim_y = data.computeSensitivity.dim_y; - dim_p = data.computeSensitivity.dim_p; +parameters = data.SWP_detection.parameters; - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +rhs = getRhsFromModelNum(datahandle, modelNum); - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +dim_y = data.computeSensitivity.dim_y; +dim_p = data.computeSensitivity.dim_p; - function_p_VDE = @(t,G) VDE_RHS_p(datahandle, sol, rhs, t, G, parameters, sensOptions); - initial_y = reshape(zeros(dim_y, dim_p), [], 1); - solVDE = integrator(function_p_VDE, tspan, initial_y, integrator_options); -end +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; + +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +function_p_VDE = @(t,G) VDE_RHS_p(datahandle, sol, rhs, t, G, parameters, sensOptions); +initial_y = reshape(zeros(dim_y, dim_p), [], 1); +solVDE = integrator(function_p_VDE, tspan, initial_y, integrator_options); +end diff --git a/toolbox/internal/sensitivities/solveVDE_Gy.m b/toolbox/internal/sensitivities/solveVDE_Gy.m index 887bf6e6..b297acc9 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gy.m +++ b/toolbox/internal/sensitivities/solveVDE_Gy.m @@ -1,25 +1,20 @@ function solVDE = solveVDE_Gy(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GY Solve the VDE for Gy in the interval tspan fixed to model modelNum and return the sol object - config = makeConfig(); - data = datahandle.getData(); +data = datahandle.getData(); - if config.removeCtrlifForSensComputation - rhs = getModelFunction(datahandle, modelNum); - else - rhs = data.integratorSettings.preprocessed_rhs; - end - parameters = data.SWP_detection.parameters; - dim_y = data.computeSensitivity.dim_y; +rhs = getRhsFromModelNum(datahandle, modelNum); +parameters = data.SWP_detection.parameters; +dim_y = data.computeSensitivity.dim_y; - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; - function_y_VDE = @(t,G) VDE_RHS_y(datahandle, sol, rhs, t, G, parameters, sensOptions); - initial_y = reshape(eye(dim_y), [], 1); - solVDE = integrator(function_y_VDE, tspan, initial_y, integrator_options); -end +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +function_y_VDE = @(t,G) VDE_RHS_y(datahandle, sol, rhs, t, G, parameters, sensOptions); +initial_y = reshape(eye(dim_y), [], 1); +solVDE = integrator(function_y_VDE, tspan, initial_y, integrator_options); +end diff --git a/toolbox/internal/solving/analyseSignature.m b/toolbox/internal/solving/analyseSignature.m index d47be841..a56097e3 100644 --- a/toolbox/internal/solving/analyseSignature.m +++ b/toolbox/internal/solving/analyseSignature.m @@ -10,7 +10,7 @@ % If ifdiff is in Filippov mode, the integration also stops upon reaching % the end of a Filippov regime. % -% If datahandle.sliding.storeSlidingInfo is true, this function stores +% If config.storeSlidingInfo is true, this function stores % detailed information about sliding mode intervals. % In particular, convex-combination parameters and the % convexified signatures are stored. @@ -39,6 +39,7 @@ % 'status': The status of the output function. % status = 0: nothing to do; status = 1: stop/halt integration +config = makeConfig(); data = datahandle.getData(); switch flag @@ -79,7 +80,7 @@ end end % sliding info - if data.sliding.storeSlidingInfo + if config.storeSlidingInfo datahandle.setData(data); convexification = analyseSignature_storeSlidingInfo(t,x,datahandle); data = datahandle.getData(); @@ -125,4 +126,4 @@ datahandle.setData(data); -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m index 43cbe795..0f2cd9d7 100644 --- a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m +++ b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m @@ -24,23 +24,10 @@ % Not in filippov mode, add NaN/entry entries. convexification.index = [convexification.index, NaN(1,length(t))]; convexification.alpha = [convexification.alpha, NaN(1,length(t))]; - convexification.signature_fplus = [convexification.signature_fplus, cell(1,length(t))]; % write empty cells - convexification.signature_fminus = [convexification.signature_fminus, cell(1,length(t))]; return end -% We don't want to redundantly store the signature many times. So we store -% it once and write a pointer in form of an index into all following entries -k = length(convexification.signature_fplus); -if k==0 || ~isscalar(convexification.signature_fplus{end}) - pointer = k+1; - convexification.signature_fplus = [convexification.signature_fplus, data.sliding.signature_fplus, repmat({pointer}, 1,length(t)-1)]; - convexification.signature_fminus = [convexification.signature_fminus, data.sliding.signature_fminus, repmat({pointer}, 1,length(t)-1)]; -else - % arrays have a pointer as the last entry, so we just add copies of it - convexification.signature_fplus = [convexification.signature_fplus, repmat(convexification.signature_fplus(end), 1,length(t))]; - convexification.signature_fminus = [convexification.signature_fminus, repmat(convexification.signature_fminus(end), 1,length(t))]; -end + % t can be vector-valued, but the Filippov-RHS can only save the alpha from % the latest evaluation, so we have to recompute it instead of reading it from @@ -54,7 +41,7 @@ % k>1 alpha = NaN(1,k); for i=1:length(t) - unused = data.sliding.filippov_rhs(datahandle, t(i), x(:,i), data.SWP_detection.parameters); %#ok + data.sliding.filippov_rhs(datahandle, t(i), x(:,i), data.SWP_detection.parameters); alpha(i) = datahandle.getData().sliding.alpha_last; end end @@ -63,4 +50,4 @@ index = datahandle.getData().sliding.index; convexification.index = [convexification.index, repmat(index, 1,length(t))]; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/chatteringGetSignatures.m b/toolbox/internal/solving/chatteringGetSignatures.m index be78a17c..359dcd8f 100644 --- a/toolbox/internal/solving/chatteringGetSignatures.m +++ b/toolbox/internal/solving/chatteringGetSignatures.m @@ -1,24 +1,21 @@ -function [signature1, signature2] = chatteringGetSignatures(datahandle, t) -% [signature1, signature2] = chatteringGetSignatures(datahandle, t) -% -% For a chattering solution, determines the signatures that appeared +function signatures = chatteringGetSignatures(datahandle, t) +% signatures = chatteringGetSignatures(datahandle, t) +% +% For a chattering solution, determines the signatures that appeared % since timepoint t. If there are more than two, throws an error since the % program is not capable of computing a Filippov solution in that case. % % % INPUT: -% datahandle: datahandle containing switching point information +% datahandle: datahandle containing switching point information % handle -% +% % t: timepoint, ignore signatures before this point % scalar % % OUTPUT: -% signature1: Signature that has been encountered. -% struct -% -% signature2: Signature that has been encountered. -% struct +% signatures: signatures that have been encountered +% 1x2 BranchingSignature data = datahandle.getData(); @@ -27,53 +24,55 @@ % Assume switching points are sorted. If not, something went very wrong... idxChatterStart = find(switchingpoints >= t, 1); -signatures = data.SWP_detection.signature; -ctrlif_indices = signatures.ctrlif_index(idxChatterStart:end); -function_indices = signatures.function_index(idxChatterStart:end); -switch_conds = signatures.switch_cond(idxChatterStart:end); - -k = length(ctrlif_indices); - -uniqueChatteringSignatures.ctrlif_index = ctrlif_indices(1); -uniqueChatteringSignatures.function_index = function_indices(1); -uniqueChatteringSignatures.switch_cond = switch_conds(1); - -for i=2:k - exists = false; - ctrlif_index = ctrlif_indices{i}; - function_index = function_indices{i}; - switch_cond = switch_conds{i}; +signatures = data.SWP_detection.signature(idxChatterStart:end); +% Make sure we don't consider any signatures from previous sliding modes. +% Sliding mode signatures are stored as 1xN BranchingSignature objects with N > 1 +isSliding = cellfun(@(x) ~isscalar(x), signatures); +if any(isSliding) + % TODO: Unsure if this case can be handled gracefully, for now abort. + throw(filippovSignatureChatteringException); +end - for j=1:length(uniqueChatteringSignatures.ctrlif_index) - ctrlif_index_ref = uniqueChatteringSignatures.ctrlif_index{j}; - function_index_ref = uniqueChatteringSignatures.function_index{j}; - switch_cond_ref = uniqueChatteringSignatures.switch_cond{j}; - if compareSignatures( ... - ctrlif_index, ctrlif_index_ref, ... - function_index, function_index_ref, ... - switch_cond, switch_cond_ref) - exists = true; +% Find unique signatures among chattering. Stop and throw if we find more than two. +uniqueChatteringSignatures = signatures(1); +for idxNew=2:length(signatures) + isNewSignature = true; + for idxOld=1:length(uniqueChatteringSignatures) + if isequal(signatures{idxNew}, uniqueChatteringSignatures{idxOld}) + % Signature already included. + isNewSignature = false; + break end end - % signature not seen yet - if ~exists - uniqueChatteringSignatures.ctrlif_index{end+1} = ctrlif_index; - uniqueChatteringSignatures.function_index{end+1} = function_index; - uniqueChatteringSignatures.switch_cond{end+1} = switch_cond; + + if isNewSignature + % New signature, add to list. + uniqueChatteringSignatures{end + 1} = signatures{idxNew}; %#ok end + + if length(uniqueChatteringSignatures) > 2 + break + end +end + + +if length(uniqueChatteringSignatures) ~= 2 + throw(invalidNumberOfChatteringSwitches); end -if length(uniqueChatteringSignatures.ctrlif_index) ~= 2 - error('IFDIFF:ChatteringInvolvesSeveralSwitches', 'Chattering does not involve exactly two signatures.\n'); +signatures = [uniqueChatteringSignatures{:}]; end -% output -signature1.ctrlif_index = uniqueChatteringSignatures.ctrlif_index{1}; -signature1.function_index = uniqueChatteringSignatures.function_index{1}; -signature1.switch_cond = uniqueChatteringSignatures.switch_cond{1}; -signature2.ctrlif_index = uniqueChatteringSignatures.ctrlif_index{2}; -signature2.function_index = uniqueChatteringSignatures.function_index{2}; -signature2.switch_cond = uniqueChatteringSignatures.switch_cond{2}; +%% Exceptions +function e = filippovSignatureChatteringException +msg = [ ... + 'Unable to determine chattering signature:\n' ... + 'One of the chattering signature candidates is a Filippov signature.\n']; +e = MException('IFDIFF:ChatteringInvolvesFilippov', msg); +end +function e = invalidNumberOfChatteringSwitches +msg = 'Chattering does not involve exactly two signatures.\n'; +e = MException('IFDIFF:ChatteringInvolvesSeveralSwitches', msg); end diff --git a/toolbox/internal/solving/compareSignatures.m b/toolbox/internal/solving/compareSignatures.m deleted file mode 100644 index 338b8205..00000000 --- a/toolbox/internal/solving/compareSignatures.m +++ /dev/null @@ -1,25 +0,0 @@ -function flag = compareSignatures(ctrlif_index1, ctrlif_index2, function_index1, function_index2, switch_cond1, switch_cond2) -% flag = compareSignatures(ctrlif_index1, ctrlif_index2, function_index1, function_index2, switch_cond1, switch_cond2) -% -% Compares two signatures. Returns true if they are equal, false if not. - -flag = true; - -if ~all(switch_cond1 == switch_cond2) || ~all(ctrlif_index1 == ctrlif_index2) - flag = false; - return -end - -if length(function_index1) ~= length(function_index2) - flag = false; - return -else - for i=1:length(function_index1) - if ~all(function_index1{i}==function_index2{i}) - flag = false; - break - end - end -end - -end \ No newline at end of file diff --git a/toolbox/internal/solving/ctrlif.m b/toolbox/internal/solving/ctrlif.m index 98e3f8d9..45e918da 100644 --- a/toolbox/internal/solving/ctrlif.m +++ b/toolbox/internal/solving/ctrlif.m @@ -1,116 +1,128 @@ function y = ctrlif(num, truepart, falsepart, ctrlif_index, function_index, datahandle) - % function for: - % - % - supervising signature during integration (.active_SWP_detection = 1) - % - forced evaluation condition during extending of the ODE solution object until switch - % - receiving the siganture at any timepoint t (.getSignature = 1) - % - forced evaluation with a given sigature for sensitivity calculation - % - % by the ctrilf one gets the condition, ctrlif_index, function_index - % those vectors are preallocated using .ctrlifCounter - - config = makeConfig(); - data = datahandle.getData(); - - condition = (num >= 0); - - switch data.caseCtrlif - case config.caseCtrlif.default - if condition - y = truepart; - else - y = falsepart; - end - - case config.caseCtrlif.forcedBranching - % active_SWP_detection: during integration used in extendODE, solveODE - % forced branching method: always return the initial true/false value - - data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; - - data.forcedBranching.switch_cond(data.forcedBranching.ctrlifCounter) = condition; - data.forcedBranching.ctrlif_index(data.forcedBranching.ctrlifCounter) = ctrlif_index; - data.forcedBranching.function_index{data.forcedBranching.ctrlifCounter} = function_index; - - - forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - - if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter - data.forcedBranching.ctrlifCounter = 0; - end - - case config.caseCtrlif.extendODEuntilSwitch - - data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; - - forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - - if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter - data.forcedBranching.ctrlifCounter = 0; - end - - case config.caseCtrlif.getSignature - data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; - - data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); - data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); - data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); - - if condition - y = truepart; - else - y = falsepart; - end - case config.caseCtrlif.getSignatureChange - data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; - - data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); - data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); - data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); - - forced_evaluation_cond = data.getSignature.switch_cond_forcedBranching(data.getSignature.ctrlifCounter); - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - case config.caseCtrlif.computeSensitivities - - % evaluate RHS with a given signature for calculating the sensitivities - - modelStage = data.computeSensitivity.modelStage; - - data.computeSensitivity.ctrlifCounter = data.computeSensitivity.ctrlifCounter + 1; - signature = data.SWP_detection.signature.switch_cond{modelStage}; - forced_evaluation_cond = signature(data.computeSensitivity.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - if length(signature) == data.computeSensitivity.ctrlifCounter - data.computeSensitivity.ctrlifCounter = 0; - end - end - - datahandle.setData(data); -end \ No newline at end of file +% function for: +% +% - supervising signature during integration (.active_SWP_detection = 1) +% - forced evaluation condition during extending of the ODE solution object until switch +% - receiving the siganture at any timepoint t (.getSignature = 1) +% - forced evaluation with a given sigature for sensitivity calculation +% +% by the ctrilf one gets the condition, ctrlif_index, function_index +% those vectors are preallocated using .ctrlifCounter + +config = makeConfig(); +data = datahandle.getData(); + +condition = (num >= 0); + +switch data.caseCtrlif + case config.caseCtrlif.default + if condition + y = truepart; + else + y = falsepart; + end + + case config.caseCtrlif.forcedBranching + % active_SWP_detection: during integration used in extendODE, solveODE + % forced branching method: always return the initial true/false value + + data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; + + data.forcedBranching.switch_cond(data.forcedBranching.ctrlifCounter) = condition; + data.forcedBranching.ctrlif_index(data.forcedBranching.ctrlifCounter) = ctrlif_index; + data.forcedBranching.function_index{data.forcedBranching.ctrlifCounter} = function_index; + + + forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + + if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter + data.forcedBranching.ctrlifCounter = 0; + end + + case config.caseCtrlif.extendODEuntilSwitch + + data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; + + forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + + if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter + data.forcedBranching.ctrlifCounter = 0; + end + + case config.caseCtrlif.getSignature + data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; + + data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); + data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); + data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); + + if condition + y = truepart; + else + y = falsepart; + end + case config.caseCtrlif.getSignatureChange + data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; + + data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); + data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); + data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); + + forced_evaluation_cond = data.getSignature.switch_cond_forcedBranching(data.getSignature.ctrlifCounter); + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + case config.caseCtrlif.computeSensitivities + + % evaluate RHS with a given signature for calculating the sensitivities + + modelStage = data.computeSensitivity.modelStage; + + data.computeSensitivity.ctrlifCounter = data.computeSensitivity.ctrlifCounter + 1; + signature = data.SWP_detection.signature{modelStage}; + + % Cannot evaluate model with sliding mode signature (i.e. multiple signatures). + if ~isscalar(signature) + throw(evaluationWithFilippovSignatureException) + end + + forced_evaluation_cond = signature.switchCond(data.computeSensitivity.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + if length(signature.switchCond) == data.computeSensitivity.ctrlifCounter + data.computeSensitivity.ctrlifCounter = 0; + end +end + +datahandle.setData(data); +end + +%% Exceptions +function e = evaluationWithFilippovSignatureException +msg = 'Cannot evaluate model in computeSensitivity mode for Filippov signature.\n'; +e = MException('IFDIFF:Ctrlif:ModelEvaluationForFilippovSignature', msg); +end diff --git a/toolbox/internal/solving/cutSteps_convexification.m b/toolbox/internal/solving/cutSteps_convexification.m index e79df427..d3f50666 100644 --- a/toolbox/internal/solving/cutSteps_convexification.m +++ b/toolbox/internal/solving/cutSteps_convexification.m @@ -16,12 +16,11 @@ function cutSteps_convexification(datahandle, t_cut) data = datahandle.getData(); convexification_t = data.sliding.convexification.t; -data.sliding.convexification.t = data.sliding.convexification.t(convexification_t<=t_cut); -data.sliding.convexification.index = data.sliding.convexification.index(convexification_t<=t_cut); -data.sliding.convexification.alpha = data.sliding.convexification.alpha(convexification_t<=t_cut); -data.sliding.convexification.signature_fplus = data.sliding.convexification.signature_fplus(convexification_t<=t_cut); -data.sliding.convexification.signature_fminus = data.sliding.convexification.signature_fminus(convexification_t<=t_cut); +mask = convexification_t <= t_cut; +data.sliding.convexification.t = data.sliding.convexification.t(mask); +data.sliding.convexification.index = data.sliding.convexification.index(mask); +data.sliding.convexification.alpha = data.sliding.convexification.alpha(mask); datahandle.setData(data); -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/evaluateRHSForSignature.m b/toolbox/internal/solving/evaluateRHSForSignature.m index fa3dc9f0..c3351162 100644 --- a/toolbox/internal/solving/evaluateRHSForSignature.m +++ b/toolbox/internal/solving/evaluateRHSForSignature.m @@ -1,5 +1,5 @@ -function dy = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index, function_index, switch_cond) -% dy = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index, function_index, switch_cond) +function dy = evaluateRHSForSignature(datahandle, t, y, p, signature) +% dy = evaluateRHSForSignature(datahandle, t, y, p, signature) % % For a given signature, evaluates the preprocessed right-hand-side in % datahandle with activated forced-branching for the given (t,y,p). @@ -16,36 +16,30 @@ % p: parameters % scalar or array % -% ctrlif_index: the signature's crtlif_index -% -% function_index: the signature's function_index -% -% switch_cond: the signature's switch_cond +% signature: signature used for forced branching evaluation +% BranchingSignature % % OUTPUT: % dy: Output of the preprocessed right-hand-side -data = datahandle.getData(); -config = makeConfig(); -% forced branching signature -data.forcedBranching.switch_cond_forcedBranching = switch_cond; -data.forcedBranching.ctrlif_index_forcedBranching = ctrlif_index; -data.forcedBranching.function_index_forcedBranching = function_index; +config = makeConfig(); +data = datahandle.getData(); +oldData = data; -% preallocate arrays to store current signature as it would be without forced branching -data.forcedBranching.switch_cond = zeros(1,length(data.forcedBranching.switch_cond_forcedBranching)); -data.forcedBranching.ctrlif_index = zeros(1,length(data.forcedBranching.ctrlif_index_forcedBranching)); -data.forcedBranching.function_index = cell(length(data.forcedBranching.function_index_forcedBranching),1); +% Set forced branching signature +data.forcedBranching.switch_cond_forcedBranching = signature.switchCond; +data.forcedBranching.ctrlif_index_forcedBranching = signature.ctrlifIndex; +data.forcedBranching.function_index_forcedBranching = signature.functionIndex; -% counter from the ctrlif +% Ensure ctrlif counter is reset data.forcedBranching.ctrlifCounter = 0; -data.caseCtrlif = config.caseCtrlif.forcedBranching; % case forced branching +% Enable forced branching mode but without storing the evaluated signature (we only care about getting output dy here). +data.caseCtrlif = config.caseCtrlif.extendODEuntilSwitch; +% Evaluate RHS with given signature datahandle.setData(data); +dy = data.integratorSettings.preprocessed_rhs(datahandle, t, y, p); -dy = datahandle.getData().integratorSettings.preprocessed_rhs(datahandle,t,y,p); - -% cleanup -data = datahandle.getData(); -data.caseCtrlif = config.caseCtrlif.default; % default case -end \ No newline at end of file +% Undo any changes made to the datahandle +datahandle.setData(oldData); +end diff --git a/toolbox/internal/solving/extendODE_filippov_regime.m b/toolbox/internal/solving/extendODE_filippov_regime.m index 7bb4e3ef..fa507b07 100644 --- a/toolbox/internal/solving/extendODE_filippov_regime.m +++ b/toolbox/internal/solving/extendODE_filippov_regime.m @@ -19,6 +19,7 @@ function extendODE_filippov_regime(datahandle) data = datahandle.getData(); data.SWP_detection.solution_until_t2 = z; data.SWP_detection.t2 = data.SWP_detection.solution_until_t2.x(end); +data.SWP_detection.x2 = data.SWP_detection.solution_until_t2.y(:, end); % write the same solution to t3 % that way, if the end of timespan is reached in filippov mode (i.e., here), @@ -26,9 +27,6 @@ function extendODE_filippov_regime(datahandle) data.SWP_detection.solution_until_t3 = z; data.SWP_detection.t3 = data.SWP_detection.solution_until_t2.x(end); -% clear filippov data, serves also as indicator for inactive filippov mode -data.sliding = extendODE_filippov_regime_cleanup(data.sliding); - datahandle.setData(data); -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/extendODEuntilSwitch.m b/toolbox/internal/solving/extendODEuntilSwitch.m index 28d96be0..c61c703d 100644 --- a/toolbox/internal/solving/extendODEuntilSwitch.m +++ b/toolbox/internal/solving/extendODEuntilSwitch.m @@ -51,6 +51,14 @@ 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); + if data.SWP_detection.t2 >= data.SWP_detection.t3 + % Cannot find a signature change between t2 and t3. + % Use t3 as the switching point, since we know that a switching event has occurred there. + data.SWP_detection.t2 = data.SWP_detection.t3; + data.SWP_detection.solution_until_t2 = data.SWP_detection.solution_until_t3; + datahandle.setData(data); + break + end datahandle.setData(data); % Check if there is a new signature. @@ -73,4 +81,4 @@ function extendODEuntilSwitch(datahandle) data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; datahandle.setData(data) -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/initDatahandleFields_Filippov.m b/toolbox/internal/solving/initDatahandleFields_Filippov.m index 3f12be04..651eb95c 100644 --- a/toolbox/internal/solving/initDatahandleFields_Filippov.m +++ b/toolbox/internal/solving/initDatahandleFields_Filippov.m @@ -11,13 +11,8 @@ sliding = extendODE_filippov_regime_cleanup(struct()); -% if true: convexification parameters and signatures are stored for every integration step -sliding.storeSlidingInfo = false; - convexification.t = []; convexification.index = []; convexification.alpha = []; -convexification.signature_fminus = {}; -convexification.signature_fplus = {}; sliding.convexification = convexification; end diff --git a/toolbox/internal/solving/initDatahandleFields_SWP_detection.m b/toolbox/internal/solving/initDatahandleFields_SWP_detection.m index 51ed2463..da529cc9 100644 --- a/toolbox/internal/solving/initDatahandleFields_SWP_detection.m +++ b/toolbox/internal/solving/initDatahandleFields_SWP_detection.m @@ -26,13 +26,10 @@ SWP_detection.function_index_t3 = []; SWP_detection.signature = {}; - SWP_detection.signature.function_index = {}; - SWP_detection.signature.ctrlif_index = {}; - SWP_detection.signature.switch_cond = {}; SWP_detection.solution_until_t1 = {}; SWP_detection.solution_until_t2 = {}; SWP_detection.solution_until_t3 = {}; SWP_detection.uniqueSwEnumeration = 0; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m b/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m index 4ce0007b..bb73940a 100644 --- a/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m +++ b/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m @@ -1,18 +1,16 @@ -function dy = slidingFilippovRHS_oneSwitch(datahandle, signature_fplus, signature_fminus, switchingFunction,t,y,p) -% Computes a Filippov right-hand-side w.r.t. the given 'ctrlifCounter', -% i.e., a convex-combination of the two models where the associated -% ctrlif evaluates to 0 and 1, respectively. +function dy = slidingFilippovRHS_oneSwitch(datahandle, signatures, switchingFunction, t, y, p) +% Computes a Filippov right-hand-side w.r.t. the given 'ctrlifCounter', +% i.e., a convex-combination of the two models where the associated +% ctrlif evaluates to 0 and 1, respectively. % Read [1] or [2, Chapter 1] for mathematical details. -% +% % INPUT: % 'datahandle' --> Datahandle for a switched system. % handle -% 'signature_fplus' --> signature corresponding to fplus in the -% sliding RHS formula -% struct -% 'signature_fminus' --> signature corresponding to fminus -% struct -% 'switchingFunction' --> Switching function of the switch that is +% 'signatures' --> signatures whose convex-combination is to be computed +% Note: first signature should have positive switching function value +% 1x2 BranchingSignature +% 'switchingFunction' --> Switching function of the switch that is % to be put in sliding mode. % function handle % 't' --> time, passed to the preprocessed RHS @@ -23,61 +21,67 @@ % scalar or array % % OUTPUT: -% 'dy' --> Value of the Filippov-RHS, determined as +% 'dy' --> Value of the Filippov-RHS, determined as % described above. % scalar or array % % Author: Michael Strik, Jun2024 % Email: michael.strik@stud.uni-heidelberg.de % michi.strik@gmail.com -% +% % [1] A.F. Filippov. Differential Equations with Discontinuous Right % Hand Side. Kluwer Academic Publishers, Dordrecht, Boston, London, 1988. -% [2] A. Meyer. Numerical Solution of Optimal Control Problems with -% Explicit and Implicit Switches. PhD thesis, Ruprecht-Karls-Universität +% [2] A. Meyer. Numerical Solution of Optimal Control Problems with +% Explicit and Implicit Switches. PhD thesis, Ruprecht-Karls-Universität % Heidelberg, 2020. -% evaluate RHS for signature where chattering switchingFunction is positive -ctrlif_index_fminus = signature_fminus.ctrlif_index; -function_index_fminus = signature_fminus.function_index; -switch_cond_fminus = signature_fminus.switch_cond; -f_minus = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index_fminus, function_index_fminus, switch_cond_fminus); - -% evaluate RHS for signature where chattering switchingFunction is negative -ctrlif_index_fplus = signature_fplus.ctrlif_index; -function_index_fplus = signature_fplus.function_index; -switch_cond_fplus = signature_fplus.switch_cond; -f_plus = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index_fplus, function_index_fplus, switch_cond_fplus); - +% Evaluate RHS for chattering signatures where switching function is positive/negative. +% Important: Assume that first signature in signatures is the positive one (has to be ensured by caller)! +f_plus = evaluateRHSForSignature(datahandle, t, y, p, signatures(1)); +f_minus = evaluateRHSForSignature(datahandle, t, y, p, signatures(2)); -FDstep = generateFDstep(length(y), length(p)); % TODO use relative finite differences (consider typical values of y and p) -n = del_f_del_y(datahandle, switchingFunction, t, y, p, FDstep.y); % n (= d/dy switchingFunction) is an outer normal of switching manifold that we want to slide on +% TODO: use relative finite differences (consider typical values of y and p) +FDstep = generateFDstep(length(y), length(p)); +n = del_f_del_y(datahandle, switchingFunction, t, y, p, FDstep.y); +n_dot_fplus = dot(n, f_plus); +n_dot_fminus = dot(n, f_minus); -n_dot_fplus = dot(n,f_plus); -n_dot_fminus = dot(n,f_minus); - -% assemble alpha -alpha = n_dot_fplus/(n_dot_fplus-n_dot_fminus); - -% SANITY CHECK -% dot(n,f_plus) and dot(n, f_minus) can't have the same sign. otherwise we -% are not in sliding mode anymore and can not rely on the previous formula for alpha +% The dot products dot(n, f_plus) and dot(n, f_minus) can't have the same sign. +% Otherwise we are not in sliding mode anymore and can not rely on the standard formula for alpha. if n_dot_fplus < 0 && n_dot_fminus < 0 - % neg: f_plus, f_minus point into {swFct<0}, continue with f_minus + % Neg: f_plus and f_minus point into {swFct<0}, continue with f_minus alpha = 1; -end - -if n_dot_fplus > 0 && n_dot_fminus > 0 - % pos: n points into {swFct>=0}, so f_plus,f_minus point into {sigma>0}. +elseif n_dot_fplus > 0 && n_dot_fminus > 0 + % Pos: f_plus and f_minus point into {swFct>0}, continue with f_plus alpha = 0; +else + % Compute alpha as usual + alpha = n_dot_fplus / (n_dot_fplus-n_dot_fminus); end -% write alpha into datahandle +% Store alpha in datahandle data = datahandle.getData(); data.sliding.alpha_last = alpha; datahandle.setData(data); -dy = alpha*f_minus + (1-alpha)*f_plus; +% Compute sliding mode RHS +if isfinite(alpha) + dy = alpha*f_minus + (1-alpha)*f_plus; +else + warnAlphaNotFinite(t, y, p); + dy = inf(size(y)); +end +end + -end \ No newline at end of file +%% Warnings +function warnAlphaNotFinite(t, y, p) +id = 'IFDIFF:Filippov:AlphaNotFinite'; +msg = [ ... + 'Filippov sliding mode RHS encountered non-finite convexification parameter at\nt=%g,\ny=[%s],\np=[%s].\n', ... + 'Returning inf value instead to force step size reduction and prevent NaN-poisoning.\n', ... + 'Consider checking the domain of your RHS for undefined regions or reducing the integrators maximum step size.' ... + ]; +warning(id, msg, t, arrayStrJoin(y, ',', '%g'), arrayStrJoin(p, ',', '%g')); +end diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m index 60ed72fe..eb9eba7c 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m @@ -58,4 +58,3 @@ function solveODE_computeSwitchingPoint(datahandle) - diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m index 1f077445..d2db0c97 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m @@ -14,8 +14,7 @@ % even number of zero crossings instead? Should probably address this % case by recursively cutting the interval and try if we can find a % subinterval [t1, t3*] with different signs on the edges. - error([ 'The switch between t=', num2str(bisection.t1), ' and t=', num2str(bisection.t3), ... - 'could not be determined. The signature did not change.\n']) + throw(noSignChangeException(bisection.t1, bisection.t3)); end @@ -23,4 +22,10 @@ % where at the edges the switching function is nonzero, % without provoking other issues? -end \ No newline at end of file +end + +%% Exceptions +function e = noSignChangeException(t1, t2) +msg = 'Unable to determine switching point: No sign change in switching function between t=%g and t=%g.\n'; +e = MException('IFDIFF:Bisection:NoSignChange', msg, t1, t2); +end diff --git a/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m b/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m index ed39fcf2..af05b5df 100644 --- a/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m +++ b/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m @@ -1,7 +1,7 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) -% Removes integration step that exceed time 't_cut' from ode sol object in +% Removes integration step that exceed time 't_cut' from ode sol object in % datahandle. -% +% % INPUT: % 'datahandle': datahandle containing the integration and switching data. % handle @@ -19,7 +19,7 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) % % Code adapted from solveODE_solution_until_t1.m -data = datahandle.getData(); +data = datahandle.getData(); % Step 1: In solution object, cut steps past time t_cut cutSteps_solution(datahandle, 'solution_until_t2', t_cut); @@ -28,16 +28,18 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) SWP_detection = data.SWP_detection; switchingpoints = cell2mat(SWP_detection.switchingpoints); -SWP_detection.switchingpoints = SWP_detection.switchingpoints(switchingpoints<=t_cut); -SWP_detection.switchingFunction = SWP_detection.switchingFunction(switchingpoints<=t_cut); -SWP_detection.jumpFunction = SWP_detection.jumpFunction(switchingpoints<=t_cut); +indicesCut = switchingpoints <= t_cut; + +SWP_detection.switchingpoints = SWP_detection.switchingpoints(indicesCut); +SWP_detection.switchingFunction = SWP_detection.switchingFunction(indicesCut); +SWP_detection.jumpFunction = SWP_detection.jumpFunction(indicesCut); t2 = data.SWP_detection.solution_until_t2.x(end); x2 = data.SWP_detection.solution_until_t2.y(:,end); SWP_detection.t2 = t2; SWP_detection.x2 = {x2, x2}; data.SWP_detection = SWP_detection; -datahandle.setData(data); +datahandle.setData(data); % adjust also signature & co. at t2 [switch_cond_t2, ctrlif_index_t2, function_index_all_t2] = ctrlif_getSignature(... datahandle, ... @@ -48,14 +50,10 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) data.SWP_detection.function_index_t2 = function_index_all_t2; % note: signature keeps length(switchingpoints)+1 entries, one entry for % the initial signature -data.SWP_detection.signature.ctrlif_index = data.SWP_detection.signature.ctrlif_index([true, switchingpoints<=t_cut]); -data.SWP_detection.signature.function_index = data.SWP_detection.signature.function_index([true, switchingpoints<=t_cut]); -data.SWP_detection.signature.switch_cond = data.SWP_detection.signature.switch_cond([true, switchingpoints<=t_cut]); +data.SWP_detection.signature = data.SWP_detection.signature([true, indicesCut]); datahandle.setData(data); % adjust convexification data cutSteps_convexification(datahandle, t_cut) - - -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/solveODE_firstTime.m b/toolbox/internal/solving/solveODE_firstTime.m index d36b6320..60dbe430 100644 --- a/toolbox/internal/solving/solveODE_firstTime.m +++ b/toolbox/internal/solving/solveODE_firstTime.m @@ -27,7 +27,7 @@ function solveODE_firstTime(datahandle) % get maxStep size by default defaultMaxStep = 0.1*abs(data.SWP_detection.tspan(1)-data.SWP_detection.tspan(2)); -data.integratorSettings.optionsForcedBranching.InitialStep = defaultMaxStep; +data.integratorSettings.optionsForcedBranching.InitialStep = defaultMaxStep; data.caseCtrlif = config.caseCtrlif.forcedBranching; datahandle.setData(data) % initialisation of nIntegrationSteps s.th. the first iteration of while is executed @@ -38,11 +38,11 @@ function solveODE_firstTime(datahandle) % while loop as long as there are only 2 timpoints or less during the integration % algorithm while nIntegrationSteps <= 2 - + % getData since the ode options may have been changed. data = datahandle.getData(); - - + + % caution: .optionsForcedBranching contains analyseSignature(...,datahandle). % switching point detection magic is done in analyseSignature(...,datahandle). z = data.integratorSettings.numericIntegrator(... @@ -50,51 +50,33 @@ function solveODE_firstTime(datahandle) data.SWP_detection.tspan,... data.SWP_detection.initialvalues, ... data.integratorSettings.optionsForcedBranching); - + data = datahandle.getData(); data.SWP_detection.solution_until_t3 = z; - + nIntegrationSteps = length(data.SWP_detection.solution_until_t3.x); - + % reduce max step size if number of timepoints not high enough if nIntegrationSteps <= 2 - + data.integratorSettings.optionsForcedBranching.InitialStep = data.integratorSettings.optionsForcedBranching.InitialStep/2; %warning('First integration step contains switch. odeextend would fail, reduce step size only for first integration step') end - + % save new MaxStep size to datahandle. - datahandle.setData(data); + datahandle.setData(data); end % get signature data.SWP_detection.t1 = data.SWP_detection.solution_until_t3.x(end-1); data.SWP_detection.t3 = data.SWP_detection.solution_until_t3.x(end); -% signature matrix -data.SWP_detection.signature.function_index{1} = data.forcedBranching.function_index_forcedBranching; -data.SWP_detection.signature.ctrlif_index{1} = data.forcedBranching.ctrlif_index_forcedBranching; -data.SWP_detection.signature.switch_cond{1} = data.forcedBranching.switch_cond_forcedBranching; -datahandle.setData(data); - +firstModelSignature = BranchingSignature( ... + data.mtreeplus{2, 1}, ... + data.forcedBranching.switch_cond_forcedBranching, ... + data.forcedBranching.ctrlif_index_forcedBranching, ... + data.forcedBranching.function_index_forcedBranching); +data.SWP_detection.signature{1} = firstModelSignature; +datahandle.setData(data); end % first integraion round done, if no switch occured, .t3 = tspan(2) - - - - - - - - - - - - - - - - - - - diff --git a/toolbox/internal/solving/solveODE_prepareNextStage.m b/toolbox/internal/solving/solveODE_prepareNextStage.m index 15416dbf..e40ded77 100644 --- a/toolbox/internal/solving/solveODE_prepareNextStage.m +++ b/toolbox/internal/solving/solveODE_prepareNextStage.m @@ -16,9 +16,8 @@ function solveODE_prepareNextStage(datahandle) % Get the new signature at t2, without forced branching [switch_cond, ctrlif_index, function_index] = ctrlif_getSignature(datahandle, t, xPlus); -data.SWP_detection.signature.switch_cond{end+1} = switch_cond; -data.SWP_detection.signature.ctrlif_index{end+1} = ctrlif_index; -data.SWP_detection.signature.function_index{end+1} = function_index; +signatureNew = BranchingSignature(data.mtreeplus{2, 1}, switch_cond, ctrlif_index, function_index); +data.SWP_detection.signature{end + 1} = signatureNew; datahandle.setData(data); end diff --git a/toolbox/internal/solving/solveODE_prepareNextStageFilippov.m b/toolbox/internal/solving/solveODE_prepareNextStageFilippov.m new file mode 100644 index 00000000..1bc3316b --- /dev/null +++ b/toolbox/internal/solving/solveODE_prepareNextStageFilippov.m @@ -0,0 +1,54 @@ +function solveODE_prepareNextStageFilippov(datahandle) +%SOLVEODE_PREPARENEXTSTAGEFILIPPOV(datahandle) +% +%Initialize datahandle fields for next round of integration after leaving Filippov sliding mode. +% +%INPUT: +% SWP_detection.t2 - Time point where the Filippov sliding mode is no longer active. +% double +% +% SWP_detection.x2 - Corresponding solution value at the exit time point. +% 1xn double +% +%OUTPUT: +% SWP_detection.switchingpoints - Filippov exit time appended as switching point. +% 1x(m+1) cell array of double +% +% SWP_detection.switchingFunction - Filippov switching function appended. +% 1x(m+1) cell array of function_handle +% +% SWP_detection.jumpFunction - Empty jump function appended (assume no jumps in Filippov). +% 1x(m+1) cell array of (function_handle | empty) +% +% SWP_detection.signature - Signature of new model at Filippov exit point appended. +% 1x(m+2) cell array of 1x? BranchingSignature +% +% sliding - Clear sliding mode fields, thus disabling Filippov sliding mode. +% struct + +data = datahandle.getData(); +swp = data.SWP_detection; +t = swp.t2; +x = swp.x2; +% First entry pre-jump, second post-jump. Here same since assume no jump function. +swp.x2 = {x, x}; + +% Exact sliding mode exit time unimportant, since transition into new model is continuous. +% Therefore, just use the last integration step (where the exit was first detected). +swp.switchingpoints{end + 1} = t; +swp.switchingFunction(end + 1) = swp.switchingFunction(end); + +% Assume no jumps in Filippov sliding mode. +swp.jumpFunction{end + 1} = []; + +% Set signature of next model. +[switch_cond, ctrlif_index, function_index] = ctrlif_getSignature(datahandle, t, x); +signatureNew = BranchingSignature(data.mtreeplus{2, 1}, switch_cond, ctrlif_index, function_index); +swp.signature{end + 1} = signatureNew; +data.SWP_detection = swp; + +% Clear Filippov data, thus disabling Filippov mode. +data.sliding = extendODE_filippov_regime_cleanup(data.sliding); + +datahandle.setData(data); +end diff --git a/toolbox/internal/solving/solveODE_setFilippovRHS.m b/toolbox/internal/solving/solveODE_setFilippovRHS.m index 0ee7526c..a518f747 100644 --- a/toolbox/internal/solving/solveODE_setFilippovRHS.m +++ b/toolbox/internal/solving/solveODE_setFilippovRHS.m @@ -1,9 +1,9 @@ function solveODE_setFilippovRHS(datahandle) -% Sets datahandle.sliding.filippovRHS to a RHS that allows to slide/run +% Sets datahandle.sliding.filippovRHS to a RHS that allows to slide/run % on the zero-manifold associated to a function that is inconsistently switching. % % INPUT: -% 'datahandle': datahandle containing the integration and +% 'datahandle': datahandle containing the integration and % switching data. % % OUTPUT: @@ -20,41 +20,36 @@ function solveODE_setFilippovRHS(datahandle) t = solveODE_backtrackChattering(datahandle); % determine the signatures involved -[signature1, signature2] = chatteringGetSignatures(datahandle, t); +signatures = chatteringGetSignatures(datahandle, t); -% only two signatures involved, so we can assume the last switching ctrlif -% to be the chattering one +% TODO: How should the index of the chattering ctrlif be determined? For now assume last ctrlif is chattering. sliding_index = datahandle.getData().SWP_detection.switchingIndices(end); -if signature1.switch_cond(sliding_index) == 1 - signature_fplus = signature1; - signature_fminus = signature2; -else - signature_fplus = signature2; - signature_fminus = signature1; +% Ensure that first signature has non-negative switching function value, i.e. switchCond is 1/true +if signatures(1).switchCond(sliding_index) ~= 1 + % Need to swap. + signatures([1 2]) = signatures([2 1]); end - % cut steps solveODE_cutSteps_solution_until_t2(datahandle, t) +data = datahandle.getData(); + % get the chattering switch's switching function: -% We ensured there was only a single switch active during chattering, +% We ensured there was only a single switch active during chattering, % so we just pick the switching function of the last switching event. -switchingFunction = datahandle.getData().SWP_detection.switchingfunctionhandles{end}; +% TODO: Is that assumption actually valid? +switchingFunction = data.SWP_detection.switchingfunctionhandles{end}; -filippov_rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch(datahandle, signature_fplus, signature_fminus, switchingFunction,t,y,p); +filippov_rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch(datahandle, signatures, switchingFunction, t, y, p); -% set filippov rhs and store some other info -data = datahandle.getData(); -data.sliding.filippov_rhs = filippov_rhs; -data.sliding.index = sliding_index; -data.sliding.signature_fplus = signature_fplus; -data.sliding.signature_fminus = signature_fminus; +% Store information about the sliding mode submodel in datahandle +data.sliding.filippov_rhs = filippov_rhs; +data.sliding.index = sliding_index; +data.SWP_detection.signature{end} = signatures; datahandle.setData(data); -% message if config.debugMode fprintf("Entered Filippov regime.\n"); end - -end \ No newline at end of file +end diff --git a/toolbox/makeConfig.m b/toolbox/makeConfig.m index d235daf1..54deab6f 100644 --- a/toolbox/makeConfig.m +++ b/toolbox/makeConfig.m @@ -30,6 +30,9 @@ % Flag to set debug mode which may display additional information and warnings among other things. config.debugMode = false; +% Flag to enable/disable storing of Filippov convexification parameters for debugging (can be costly!). +config.storeSlidingInfo = false; + % Flag to enable/disable reusing of switching functions that already exist based on the signature of a switch. % Always keep this enabled for best performance. Can be useful to disable for debugging and tests. config.reuseSwitchingFunctions = true; diff --git a/toolbox/solveODE.m b/toolbox/solveODE.m index bfa952cc..0212e33c 100644 --- a/toolbox/solveODE.m +++ b/toolbox/solveODE.m @@ -60,14 +60,21 @@ if checkForSwitchingIndices(datahandle) error('IFDIFF:switchInFilippov', 'Switching during or right after Filippov mode.'); end - if datahandle.getData().SWP_detection.t2 == tspan(2) - break; % reached end of time span - else - if config.debugMode - fprintf('Left Filippov regime.\n'); - end - solveODE_prepareNextStage(datahandle); + + % Filippov sliding mode integration has stopped + data = datahandle.getData(); + if data.SWP_detection.t2 == tspan(2) + % Reached end of integration timespan. + data.sliding = extendODE_filippov_regime_cleanup(data.sliding); + datahandle.setData(data); + break; + end + + % Left the Filippov sliding mode. Must continue integration with new model. + if config.debugMode + fprintf('Left Filippov regime.\n'); end + solveODE_prepareNextStageFilippov(datahandle); end % extend solution object from t2 ongoing until the next switch occurs