Help with code for impulsive optimal control problem

Question

I am using Mathematica 13.3 to numerically solve/generate numerical plots for my impulsive control problem.

The optimal control problem is:

T is time from0 to T where T is the terminal time

K(t), B(t) and M(t) are state variables

w(t) is a control variable

a(ti) is the impulsive control variable at time ti,

a(ti) \in [0,1] for i=1,2,…,N

ggamma, ttheta, ddelta1, ddelta 2, c1 and c2 are constants

K(T)=M(T)=0 B(T)>0

K'(t)=ggamma*K(t)w(t)-tthetaM(t)

B’(t)=ggammaK(t)+tthetaM(t)*B(t)

M’(t)=M(t)-ggamma*w(t)

M(ti)=M(ti-)+a(ti)M(ti)*ddelta1

K(ti)=K(ti-)-a(ti)K(ti)*ddelta2

Objective: maximize B(T)-integral from 0 to T of c1*(w(t))^2dt-sum i=1 to N of c2*a(ti)

 (* Define the constants *)
gamma = 1.0;'
theta = 2.0;'
delta1 = 0.1;'
delta2 = 0.2;'
c1 = 0.5;'
c2 = 0.3;'
T = 5.0; (* Terminal time *)'

(* Define the impulsive changes in M[t] *)
impulseChanges[t_] := Module[{tiValues, impValues, result},
   tiValues = {1.0, 2.0, 3.0}; (* Example impulsive time instants *)
   impValues = {0.2, 0.1, 0.3}; (* Corresponding impulsive control values *)
   result = 0;
   Do[
    If[t == tiValues[[i]],
     result = result + impValues[[i]] M[tiValues[[i]]] (delta1 - delta2)],
    {i, Length[tiValues]}];
   result
   ];

(* Define the system of differential equations *)
diffeqs = {
   K'[t] == gamma K[t] w[t] - theta M[t],
   B'[t] == gamma K[t] + theta M[t] B[t],
   M'[t] == M[t] - gamma w[t]
   };

(* Define the impulsive controls *)
impulseControls = {1.0, 0.5, 0.8}; (* Example impulsive control values *)

(* Define the initial values and derivatives *)
initialValues = {K[0], B[0], M[0]} == {10, 20, 15};
initialDerivatives = {K'[0], B'[0], M'[0]} == {0, 0, 0}; (* Adjust as needed *)

(* Define the objective function to be maximized *)
objective = 
  B[T] - Integrate[c1 w[t]^2, {t, 0, T}] - 
   Sum[c2 impulseControls[[i]] /. t -> ti, {i, Length[impulseControls]}];

(* Solve the system of differential equations numerically using the "ImplicitRungeKutta" method *)
sol = NDSolveValue[
   Join[diffeqs, initialValues, initialDerivatives],
   {K, B, M},
   {t, 0, T},
   Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 5},
   MaxSteps -> Infinity
   ];

(* Find the optimal control trajectory w[t] using optimization *)
wOptimal[t_?NumericQ] := 
  w[t] /. NMaximize[objective, w[t]][[2]];

(* Evaluate the optimal control and state trajectories *)
optimalControls = Table[wOptimal[t], {t, 0.0, T, 0.1}];
stateTrajectories = {K[t], B[t], M[t]} /. sol;

optimalControls
stateTrajectories

I am getting many error messages, like these ones:

NDSolveValue::deqn: Equation or list of equations expected instead of Join[{0==-2. M[t]+1. K[t] w[t],0==1. K[t]+2. B[t] M[t],0==M[t]-1. w[t]},{K[0],B[0],M[0]}=={10,20,15},True] in the first argument Join[{0==-2. M[t]+1. K[t] w[t],0==1. K[t]+2. B[t] M[t],0==M[t]-1. w[t]},{K[0],B[0],M[0]}=={10,20,15},True].

Integrate::ilim: Invalid integration variable or limit(s) in {0.,0,5.}.

General::stop: Further output of Integrate::ilim will be suppressed during this calculation.

NMaximize::nnum: The function value 0.69 -B[5.]+!(*SubsuperscriptBox[([Integral]), (0), (5.\)]\(0.34366403146382957 [DifferentialD]0)) is not a number at {w[0.]} = {-0.829053}.

ReplaceAll::reps: {w[0.]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

ReplaceAll::reps: {w[0.1]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

I would appreciate any help with fixing my code!

Thank you very much!