How to include impulses in differential equations?

Question

I have two coupled differential equations as follows:

$$S'(t) = - \frac{a}{1 + B(t)} S(t),$$

$$B'(t) = \frac{c}{1 + S(t)} B(t) (1 - B(t)) - d B^2(t) \Big( \frac{1 - B(t)}{B(t)} \Big)^n,$$

where $a$, $c$, $d$, and $n$ are constants to be determined from the data. And the datasets are:

Sdata = {{0, 9.74},{3, 4.92},{6, 8.29},{9, 5.54},{15, 2.08},{18, 1.38},{21, 1.99},{24, 0.893}};
Bdata = {{0, 0.915094},{3, 0.736097},{6, 0.793694},{9, 0.833664},{15, 1},{18, 0.99578},{21, 0.897964},{24, 0.214499}};

We fit the differential equations to the datasets and obtain:

Clear[fitness]
fitness[a00_?NumberQ, c00_?NumberQ, d00_?NumberQ, n00_?NumberQ, s0_?NumberQ, b0_?NumberQ] := 
Module[{ss, bb, parms, odes, ODES, solode, cinits, ssbb, lambda = 4.75, t}, 
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, bb'[t] - c00  bb[t] (1 - bb[t])/(1 + ss[t]) + d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == 0};
cinits = {ss[0] == s0, bb[0] == b0};
ODES = Join[odes, cinits];
solode = Quiet@NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
Return[
Sum[Norm[{y1[[k]], lambda y2[[k]]} - 
Evaluate[{ss[t], lambda bb[t]} /. solode /. {t -> tk[[k]]}]], {k, 1, 8}]]]
tk = {0, 3, 6, 9, 15, 18, 21, 24};
y1 = {9.74, 4.92, 8.29, 5.54, 2.08, 1.38, 1.99, 0.893};
y2 = {0.915094, 0.736097, 0.793694, 0.833664, 1, 0.99578, 0.897964, 0.214499};
restrs = {0 < a00 < 0.9 && 0 < c00 < 0.9 && 0.1 < d00 < 0.8 && 1.15 < n00 < 1.2 && 8. < s0 < 10 && 0.01 < b0 < 0.99999999};
vars = {a00, c00, d00, n00, s0, b0};
sol = Quiet@
FindMinimum[
Join[{fitness[a00, c00, d00, n00, s0, b0]}, 
restrs], {{a00, 0.13}, {c00, 0.13}, {d00, 0.24}, {n00, 0.25}, {s0, 9}, {b0, 1}}, 
StepMonitor :> Print[{a00, c00, d00, n00, s0, b0, 
fitness[a00, c00, d00, n00, s0, b0]}], MaxIterations -> 10, 
WorkingPrecision -> 30]
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, bb'[t] - c00  bb[t] (1 - bb[t])/(1 + 
ss[t]) + d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == 0} /. sol[[2]];
cinits = {ss[0] == s0, bb[0] == b0} /. sol[[2]];
ODES = Join[odes, cinits];
solode = NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
gr0 = Plot[Evaluate[ss[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr10 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Red, 
PlotLabel -> "S(t)"];
gr2 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue];
Show[gr10, gr0]
gr00 = Plot[Evaluate[bb[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr02 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue, 
PlotLabel -> "B(t)"];
Show[gr02, gr00]

I'm much more interested to capture the behavior of $B(t)$ more accurately. The behavior of $B$ data (blue points) is: First it decreases, to have a minimum at $t = 3$, then increases and has a maximum at $t = 15$, and then suddenly it falls. The fitted curve is more or less flat till $t = 15$. It seems to me that coupled two-components first-order differential equations cannot capture the behavior of $B$ data, and one needs either higher derivatives or more components (at least three). However, I don't have any other information about my system to be included in the ODE's. I searched around a little and found the concept of impulsive differential equations with one example here with the associated answer here, which it seems it can be used also in my case. The main idea is to put Dirac deltas and force jumps in the solutions of differential equations, which for my case it means around $t = 3$ and $t = 15$.

To this end, we plug in -DiracDelta[t - 3] + DiracDelta[t - 15] in the right hand side of the differential equation for $B(t)$, two places in the above code, to do fitting again. However, Mathematica runs forever on my PC (4 GB RAM). I'd appreciate it if anyone could see if they can produce results on their PC.

EDIT 1

As Ulrich Neumann suggested in the comments, we approximate Dirac delta as dirac = Function[{t, eps}, Exp[-(t^2/(2 eps))]/Sqrt[2 Pi eps]], and plug in the term: -u3 dirac[t - 3, .1] + u15 dirac[t - 15, .1], with two new parameters $u_3$ and $u_{15}$, in the right-hand side of the differential equation for $B(t)$, thus we have:

Clear[fitness]
Sdata = {{0, 9.74}, {3, 4.92}, {6, 8.29}, {9, 5.54}, {15, 2.08},{18, 1.38}, {21, 1.99}, {24, 0.893}};
Bdata = {{0, 0.915094}, {3, 0.736097}, {6, 0.793694}, {9, 0.833664}, {15, 1}, {18, 0.99578}, {21, 0.897964}, {24, 0.214499}};
dirac = Function[{t, eps}, Exp[-(t^2/(2 eps))]/Sqrt[2 Pi eps]];
fitness[a00_?NumberQ, c00_?NumberQ, d00_?NumberQ, n00_?NumberQ, 
s0_?NumberQ, b0_?NumberQ, u3_?NumberQ, u15_?NumberQ] := 
Module[{ss, bb, parms, odes, ODES, solode, cinits, ssbb,lambda = 4.75, t}, 
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, 
bb'[t] - c00 bb[t] (1 - bb[t])/(1 + ss[t]) + d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == -u3 dirac[t - 3, .1] + 
u15 dirac[t - 15, .1]};
cinits = {ss[0] == s0, bb[0] == b0};
ODES = Join[odes, cinits];
solode = Quiet@NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
Return[
Sum[Norm[{y1[[k]], lambda y2[[k]]} - Evaluate[{ss[t], lambda bb[t]} /. 
solode /. {t -> tk[[k]]}]], {k, 1, 8}]]]
tk = {0, 3, 6, 9, 15, 18, 21, 24};
y1 = {9.74, 4.92, 8.29, 5.54, 2.08, 1.38, 1.99, 0.893};
y2 = {0.915094, 0.736097, 0.793694, 0.833664, 1, 0.99578, 0.897964, 0.214499};
restrs = {0 < a00 < 0.9 && 0 < c00 < 0.9 && 0.1 < d00 < 0.8 && 
1.15 < n00 < 1.2 && 8. < s0 < 10 && 0.01 < b0 < 0.99999999 && 
0.1 < u3 < 0.15 && 0.2 < u15 < 0.3};
vars = {a00, c00, d00, n00, s0, b0, u3, u15};
sol = Quiet@FindMinimum[Join[{fitness[a00, c00, d00, n00, s0, b0, 
u3, u15]},restrs], {{a00, 0.13}, {c00, 0.13}, {d00, 0.24}, {n00, 0.25}, {s0, 9}, {b0, 1}, {u3, 0.1}, {u15, 0.1}}, 
StepMonitor :> 
Print[{a00, c00, d00, n00, s0, b0, u3, u15, fitness[a00, c00, d00, n00, s0, b0, u3, u15]}], MaxIterations -> 10, WorkingPrecision -> 30]
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, 
bb'[t] - c00 bb[t] (1 - bb[t])/(1 + ss[t]) + d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == -u3 dirac[t - 3, .1] + 
u15 dirac[t - 15, .1]} /. sol[[2]];
cinits = {ss[0] == s0, bb[0] == b0} /. sol[[2]];
ODES = Join[odes, cinits];
solode = NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
gr0 = Plot[Evaluate[ss[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr10 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Red, 
PlotLabel -> "S(t)"];
gr2 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue];
Show[gr10, gr0]
gr00 = Plot[Evaluate[bb[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr02 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue, 
PlotLabel -> "B(t)"];
Show[gr02, gr00]

So, it seems Dirac delta is not the best approach for modeling impulsive differential equations.

Any idea for forcing impulses is really appreciated!

EDIT 2

With only one Dirac delta at $t = 3$, the fitting is a little better, but still cannot capture the behavior of the data the end:

Clear[fitness]
Sdata = {{0, 9.74}, {3, 4.92}, {6, 8.29}, {9, 5.54}, {15, 2.08},{18, 1.38}, {21, 1.99}, {24, 0.893}};
Bdata = {{0, 0.915094}, {3, 0.736097}, {6, 0.793694}, {9, 0.833664}, {15, 1}, {18, 0.99578}, {21, 0.897964}, {24, 0.214499}};
dirac = Function[{t, eps}, Exp[-(t^2/(2 eps))]/Sqrt[2 Pi eps]];
fitness[a00_?NumberQ, c00_?NumberQ, d00_?NumberQ, n00_?NumberQ, s0_?NumberQ, b0_?NumberQ, u3_?NumberQ, u15_?NumberQ] := 
Module[{ss, bb, parms, odes, ODES, solode, cinits, ssbb, lambda = 4.75, t}, 
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, 
bb'[t] - c00 bb[t] (1 - bb[t])/(1 + ss[t]) - d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == -u3 dirac[t - 3, .5]};
cinits = {ss[0] == s0, bb[0] == b0};
ODES = Join[odes, cinits];
solode = Quiet@NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
Return[Sum[Norm[{y1[[k]], lambda y2[[k]]} - Evaluate[{ss[t], lambda bb[t]} /. solode /. {t -> tk[[k]]}]], {k, 1, 8}]]]
tk = {0, 3, 6, 9, 15, 18, 21, 24};
y1 = {9.74, 4.92, 8.29, 5.54, 2.08, 1.38, 1.99, 0.893};
y2 = {0.915094, 0.736097, 0.793694, 0.833664, 1, 0.99578, 0.897964, 0.214499};
restrs = {0 < a00 < 0.9 && 0 < c00 < 0.9 && 0.1 < d00 < 0.8 && 
1.15 < n00 < 1.2 && 8. < s0 < 10 && 0.01 < b0 < 0.99999999 && 
0.1 < u3 < 0.15 && 0.2 < u15 < 0.3};
vars = {a00, c00, d00, n00, s0, b0, u3, u15};
sol = Quiet@FindMinimum[
Join[{fitness[a00, c00, d00, n00, s0, b0, u3, u15]}, 
restrs], {{a00, 0.13}, {c00, 0.13}, {d00, 0.24}, {n00, 0.25}, {s0, 9}, {b0, 1}, {u3, 0.1}, {u15, 0.1}}, 
StepMonitor :> 
Print[{a00, c00, d00, n00, s0, b0, u3, u15, 
fitness[a00, c00, d00, n00, s0, b0, u3, u15]}], 
MaxIterations -> 10, WorkingPrecision -> 30]
odes = {ss'[t] + a00/(1 + bb[t]) ss[t] == 0, 
bb'[t] - c00 bb[t] (1 - bb[t])/(1 + ss[t]) - d00 bb[t]^2 ((1 - bb[t])/bb[t])^n00 == -u3 dirac[t - 3, .5]} /. sol[[2]];
cinits = {ss[0] == s0, bb[0] == b0} /. sol[[2]];
ODES = Join[odes, cinits];
solode = NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
gr0 = Plot[Evaluate[ss[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr10 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Red, 
PlotLabel -> "S(t)"];
gr2 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue];
Show[gr10, gr0]
gr00 = Plot[Evaluate[bb[t] /. solode], {t, tk[[1]], tk[[8]]}];
gr02 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue, 
PlotLabel -> "B(t)"];
Show[gr02, gr00]