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/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..0531bbdb --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m @@ -0,0 +1,213 @@ +%% Setup +% Time span +tspan = [0, 100]; + +% Integrator +intIfdiff = @ode45; +intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12, 'MaxStep', 0.5); +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 +data = datahandle.getData(); +data.sliding.storeSlidingInfo = true; +datahandle.setData(data); + +solIfdiff = solveODE(datahandle, tspan, x0, p); + + +%% 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..781be0e7 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,34 @@ this.functionIndex = val; this = this.updateCache(); end + + %% Comparison operators + % Equality operator overloading + function tf = eq(obj1, obj2) + tf = false; + % Objects should be of same class, no conversions. + if ~strcmp(class(obj1), class(obj2)) + return + end + + % Check each property + if ~strcmp(obj1.rhsName, obj2.rhsName) + return + end + if ~isequal(obj1.switchCond, obj2.switchCond) + return + end + if ~isequal(obj1.ctrlifIndex, obj2.ctrlifIndex) + return + end + if ~isequal(obj1.functionIndex, obj2.functionIndex) + return + end + tf = true; + end + + function tf = ne(obj1, obj2) + tf = ~eq(obj1, obj2); + end end end diff --git a/toolbox/internal/rhsPreprocessing/createDatahandle.m b/toolbox/internal/rhsPreprocessing/createDatahandle.m index 143fe1a1..daa2f331 100644 --- a/toolbox/internal/rhsPreprocessing/createDatahandle.m +++ b/toolbox/internal/rhsPreprocessing/createDatahandle.m @@ -58,5 +58,7 @@ data.caseCtrlif = config.caseCtrlif.default; +data.sliding.storeSlidingInfo = config.storeSlidingInfo; + datahandle = makeClosure(data); 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..eda9de82 --- /dev/null +++ b/toolbox/internal/sensitivities/getRhsFromModelNum.m @@ -0,0 +1,48 @@ +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 iscell(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{1}, ... + signature{2}, ... + 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_storeSlidingInfo.m b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m index 43cbe795..7c301d6b 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 @@ -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..709d74a1 100644 --- a/toolbox/internal/solving/chatteringGetSignatures.m +++ b/toolbox/internal/solving/chatteringGetSignatures.m @@ -1,22 +1,22 @@ function [signature1, signature2] = chatteringGetSignatures(datahandle, t) % [signature1, signature2] = chatteringGetSignatures(datahandle, t) -% -% For a chattering solution, determines the signatures that appeared +% +% 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 @@ -27,53 +27,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 Filippov modes. +% Filippov signatures are stored as cells with two signatures. +isFilippov = cellfun(@iscell, signatures); +if any(isFilippov) + % 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 signatures{idxNew} == signatures{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.ctrlif_index) ~= 2 - error('IFDIFF:ChatteringInvolvesSeveralSwitches', 'Chattering does not involve exactly two signatures.\n'); + +if length(uniqueChatteringSignatures) ~= 2 + throw(invalidNumberOfChatteringSwitches); end % output -signature1.ctrlif_index = uniqueChatteringSignatures.ctrlif_index{1}; -signature1.function_index = uniqueChatteringSignatures.function_index{1}; -signature1.switch_cond = uniqueChatteringSignatures.switch_cond{1}; +[signature1, signature2] = uniqueChatteringSignatures{:}; +end -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..7d36a81c 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 Filippov signature (i.e. two signatures). + if iscell(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..d96ec6c2 100644 --- a/toolbox/internal/solving/evaluateRHSForSignature.m +++ b/toolbox/internal/solving/evaluateRHSForSignature.m @@ -39,6 +39,7 @@ % counter from the ctrlif data.forcedBranching.ctrlifCounter = 0; +oldCaseCtrlif = data.caseCtrlif; data.caseCtrlif = config.caseCtrlif.forcedBranching; % case forced branching datahandle.setData(data); @@ -47,5 +48,6 @@ % cleanup data = datahandle.getData(); -data.caseCtrlif = config.caseCtrlif.default; % default case -end \ No newline at end of file +data.caseCtrlif = oldCaseCtrlif; +datahandle.setData(data); +end diff --git a/toolbox/internal/solving/extendODE_filippov_regime.m b/toolbox/internal/solving/extendODE_filippov_regime.m index 7bb4e3ef..2e82f52c 100644 --- a/toolbox/internal/solving/extendODE_filippov_regime.m +++ b/toolbox/internal/solving/extendODE_filippov_regime.m @@ -26,9 +26,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..9b626a96 100644 --- a/toolbox/internal/solving/extendODEuntilSwitch.m +++ b/toolbox/internal/solving/extendODEuntilSwitch.m @@ -51,6 +51,11 @@ 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 + % Okay, just use t3. + data.SWP_detection.t2 = data.SWP_detection.t3; + break + end datahandle.setData(data); % Check if there is a new signature. @@ -73,4 +78,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/extendODEuntilSwitch_Filippov.m b/toolbox/internal/solving/extendODEuntilSwitch_Filippov.m new file mode 100644 index 00000000..abff8e75 --- /dev/null +++ b/toolbox/internal/solving/extendODEuntilSwitch_Filippov.m @@ -0,0 +1,99 @@ +function extendODEuntilSwitch_Filippov(datahandle) +%EXTENDODEUNTILSWITCH Extend ODE solution until switch and set new signature. +% EXTENDODEUNTILSWITCH(datahandle) +% +% Extend the ODE solution from the previous time point (t1) until the switching time (t2) using forced branching. +% Note that t2 has to be part of the new model (i.e. have a new signature). To achieve this, the numerically +% computed switching time (t2), which could still be part of the old model due to inaccuracies, +% may be increased by this function. +% +% INPUT +% Accepts input via datahandle providing the fields +% - SWP_detection.switch_cond_t1/t2 : switch condition before and after the numerically computed switch +% - SWP_detection.t2 : numerically computed switching time +% - SWP_detection.solution_until_t1 : ODE solution before the switch +% +% OUTPUT +% Writes results back to datahandle in the fields +% - SWP_detection.t2 : updated switching time (guaranteed to be part of a new model) +% - SWP_detection.solution_until_t2 : ODE solution until the updated switching time +% - SWP_detection.switchingpoints : new switching time appended +% - SWP_detection.signature : new signature appended + + +config = makeConfig(); +data = datahandle.getData(); + +% Check if the numerically computed switching point is part of the new model. +extendODEuntilSwitch_t1_to_t2(datahandle); +done = isFilippovExit(datahandle); + +% The numerically computed switching point may still be part of the old model due to inaccuracies. +% In that case, slightly increment the switching point and reintegrate with forced branching. +% Repeat this process until a new signature is detected. +t2FromRootFinding = data.SWP_detection.t2; +baseOffset = 16*eps(data.SWP_detection.t2); +iter = 0; +while ~done + data = datahandle.getData(); + + % Throw error if error of switching point obtained from root finding relative to tspan exceeds threshold. + % Note: No risk of division by zero since we need to make at least one integration step to detect a switch. + switchingPointError = abs((data.SWP_detection.t2 - t2FromRootFinding) ... + / (data.SWP_detection.tspan(end) - data.SWP_detection.tspan(1))); + if switchingPointError > config.switchingPointErrorThreshold + errorMsg = 'Relative error of numerically computed switching point exceeds threshold: %.16g > %.16g\n'; + error('IFDIFF:switchingPointErrorThreshold', errorMsg, switchingPointError, config.switchingPointErrorThreshold); + end + + % Increase t2 and integrate again starting from t1. + % 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 + % Okay, just use t3. + data.SWP_detection.t2 = data.SWP_detection.t3; + break + end + datahandle.setData(data); + + % Check if there is a new signature. + extendODEuntilSwitch_t1_to_t2(datahandle); + done = isFilippovExit(datahandle); + + iter = iter + 1; +end + +data = datahandle.getData(); + +% Log events related to searching for the switch. +if config.debugMode + switchingPointInfo = "DEBUG: Detected switch after %d iterations at time point t=%.16g\n"; + fprintf(switchingPointInfo, iter, data.SWP_detection.t2); +end + +% Add new switching point to history. +data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; + +datahandle.setData(data) +end + +%% Helpers +function updateAlpha(datahandle) +data = datahandle.getData(); +t = data.SWP_detection.t2; +x = deval(data.SWP_detection.solution_until_t2, t); +unused = data.sliding.filippov_rhs(datahandle, t, x, data.SWP_detection.parameters); %#ok +end + +function tf = isFilippovExit(datahandle) +updateAlpha(datahandle); + +data = datahandle.getData(); +[abstol, reltol] = getIntegratorTolerances(data.integratorSettings.options); +state = data.SWP_detection.solution_until_t2.x(end); +alpha_tol = min(abstol, reltol*norm(state)); +alpha = data.sliding.alpha_last; + +tf = alpha <= alpha_tol || alpha >= 1-alpha_tol; +end diff --git a/toolbox/internal/solving/initDatahandleFields.m b/toolbox/internal/solving/initDatahandleFields.m index 344e95ae..802aedb8 100644 --- a/toolbox/internal/solving/initDatahandleFields.m +++ b/toolbox/internal/solving/initDatahandleFields.m @@ -16,7 +16,7 @@ function initDatahandleFields(datahandle, tspan, initialvalues, parameters) data.forcedBranching = initDatahandleFields_forcedBranching(); data.getSignature = initDatahandleFields_getSignature(); data.computeSensitivity = initDatahandleFields_computeSensitivity(); -data.sliding = initDatahandleFields_Filippov(); +data.sliding = initDatahandleFields_Filippov(data.sliding); datahandle.setData(data); diff --git a/toolbox/internal/solving/initDatahandleFields_Filippov.m b/toolbox/internal/solving/initDatahandleFields_Filippov.m index 3f12be04..7ab85e3f 100644 --- a/toolbox/internal/solving/initDatahandleFields_Filippov.m +++ b/toolbox/internal/solving/initDatahandleFields_Filippov.m @@ -1,23 +1,18 @@ -function sliding = initDatahandleFields_Filippov() +function sliding = initDatahandleFields_Filippov(sliding) %sliding = INITDATAHANDLEFIELDS_FILIPPOV() % %Initialize struct for Filippov mode handling in IFDIFF. % -%OUTPUT: +%INPUT/OUTPUT: % sliding - Struct used for Filippov mode handling in IFDIFF. % struct % %See also INITDATAHANDLEFIELDS. -sliding = extendODE_filippov_regime_cleanup(struct()); - -% if true: convexification parameters and signatures are stored for every integration step -sliding.storeSlidingInfo = false; +sliding = extendODE_filippov_regime_cleanup(sliding); 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..7b95a7fc 100644 --- a/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m +++ b/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m @@ -38,15 +38,15 @@ % 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; +ctrlif_index_fminus = signature_fminus.ctrlifIndex; +function_index_fminus = signature_fminus.functionIndex; +switch_cond_fminus = signature_fminus.switchCond; 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; +ctrlif_index_fplus = signature_fplus.ctrlifIndex; +function_index_fplus = signature_fplus.functionIndex; +switch_cond_fplus = signature_fplus.switchCond; f_plus = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index_fplus, function_index_fplus, switch_cond_fplus); @@ -80,4 +80,4 @@ dy = alpha*f_minus + (1-alpha)*f_plus; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m index 60ed72fe..e9faa458 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m @@ -8,7 +8,7 @@ function solveODE_computeSwitchingPoint(datahandle) % preallocate for switchingIntervals -numberofswitchingfunctions = length(data.SWP_detection.switchingIndices); +numberofswitchingfunctions = length(data.SWP_detection.switchingfunctionhandles); switchingpoints = zeros(1, numberofswitchingfunctions); data.caseCtrlif = config.caseCtrlif.default; % find zero crossing of the switching function in the time interval diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m index 1f077445..9e4dc11f 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m @@ -2,25 +2,57 @@ % Function that ensures that the following bisection algorithm can run % without issues. % In particular, it ensure that bisection.switchingFunction doesn't -% evaluate to zero at bisection.t1 and that it changes sign between +% evaluate to zero at bisection.t1 and that it changes sign between % bisection.t1 and bisection.t3. % % FUNCTION NEEDS REVISION. +% Sign change, everything is fine. +if sign(bisection.sw1) ~= sign(bisection.sw3) + return +end -if sign(bisection.sw1) == sign(bisection.sw3) - % no switch has occurred +tolZero = 1e-10; +% Both interval points are far from zero, so inputs are bad, nothing we can do. +if abs(bisection.sw1) >= tolZero && abs(bisection.sw3) >= tolZero % TODO: Can we be certain that the switching function doesn't have an % 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 +% At least one of the interval points is close to zero. +% Might be able to create a zero crossing via small shift in switching function. +shift = 2 * min(abs(bisection.sw1), abs(bisection.sw3)); +% If points positive shift down, else shift up. +if bisection.sw1 >= 0 + shift = -shift; +end +bisection.sw1 = bisection.sw1 + shift; +bisection.sw3 = bisection.sw3 + shift; +% Check sign change again after shift +if sign(bisection.sw1) == sign(bisection.sw3) + throw(noSignChangeException(bisection.t1, bisection.t3)); +end +bisection.switchingFunction = @(dh, t, x, p) bisection.switchingFunction(dh, t, x, p) + shift; +warnApplyingShift; -% Open question: What if bisection.sw1 == 0? Can we find an interval +% Open question: What if bisection.sw1 == 0? Can we find an interval % where at the edges the switching function is nonzero, % without provoking other issues? +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 -end \ No newline at end of file +%% Warnings +function warnApplyingShift +msg = [ + 'Interval points for switching point computation have same sign, but close to zero.\n', ... + 'Applying small shift to switching function in an attempt to restore sign change.\n']; +warning('IFDIFF:Bisection:ApplyShift', msg); +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_setFilippovRHS.m b/toolbox/internal/solving/solveODE_setFilippovRHS.m index 0ee7526c..3cf1012a 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: @@ -25,7 +25,7 @@ function solveODE_setFilippovRHS(datahandle) % only two signatures involved, so we can assume the last switching ctrlif % to be the chattering one sliding_index = datahandle.getData().SWP_detection.switchingIndices(end); -if signature1.switch_cond(sliding_index) == 1 +if signature1.switchCond(sliding_index) == 1 signature_fplus = signature1; signature_fminus = signature2; else @@ -33,28 +33,28 @@ function solveODE_setFilippovRHS(datahandle) signature_fminus = signature1; end - % cut steps solveODE_cutSteps_solution_until_t2(datahandle, t) +data = datahandle.getData(); + +% Store signature of Filippov model as combined signatures for sensitivities. +data.SWP_detection.signature{end} = {signature_fplus, signature_fminus}; + % 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}; +switchingFunction = data.SWP_detection.switchingfunctionhandles{end}; filippov_rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch(datahandle, signature_fplus, signature_fminus, 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; + 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..0f653d33 100644 --- a/toolbox/solveODE.m +++ b/toolbox/solveODE.m @@ -60,14 +60,41 @@ if checkForSwitchingIndices(datahandle) error('IFDIFF:switchInFilippov', 'Switching during or right after Filippov mode.'); end - if datahandle.getData().SWP_detection.t2 == tspan(2) + + data = datahandle.getData(); + if data.SWP_detection.t2 == tspan(2) + % clear filippov data, serves also as indicator for inactive filippov mode + data.sliding = extendODE_filippov_regime_cleanup(data.sliding); + datahandle.setData(data); break; % reached end of time span - else - if config.debugMode - fprintf('Left Filippov regime.\n'); - end - solveODE_prepareNextStage(datahandle); end + + % Left the Filippov model and continue integration with new model: + % Need to set switching point, switching function and new signature. + if config.debugMode + fprintf('Left Filippov regime.\n'); + end + + solveODE_solution_until_t1(datahandle) + + data = datahandle.getData(); + % Switching function for entry into Filippov regime. + data.SWP_detection.switchingfunctionhandles = data.SWP_detection.switchingFunction(end); + data.SWP_detection.switchingIndices = 1; % Doesn't matter, but required for computeSwitchingPoint + % Assume no jumps in Filippov regime. + data.SWP_detection.jumpFunction{end + 1} = []; + datahandle.setData(data); + solveODE_computeSwitchingPoint(datahandle); + + % Compute actual switching point by extending solution until alpha changes. + extendODEuntilSwitch_Filippov(datahandle); + + % clear filippov data, serves also as indicator for inactive filippov mode + data = datahandle.getData(); + data.sliding = extendODE_filippov_regime_cleanup(data.sliding); + datahandle.setData(data); + + solveODE_prepareNextStage(datahandle); end % extend solution object from t2 ongoing until the next switch occurs