How can I solve following cooperative differential game? Consider a game-theoretic model of pollution control. There are 2 players join in the game, N = {1, 2}. Each player has an industrial production site. It is assumed
that the production is proportional to the pollutions $u_i$. Thus, the strategy of a player is to choose the amount of pollutions emited to the atmosphere, $u_i \in [0; b_i]$. In this example the solution will be
considered in the class of open-loop strategies $u_i(t)$.
The dynamics of the total amount of pollution x(t) is described by
$\dot x=u_1+u_2, x(t_0)=x_0.$
The payoff of the i-th player is defined as
$\int_{0}^{T} (b_i −\frac{1}{2}u_i)u_i − d_ix)dt, i = 1, 2.$
We assume that the terminal cost is zero.
Each player is to maximize their total cost, then the optimization problem is as follows:
$\sum_{i=1}^{2}\int_{0}^{T}(b_i-\frac{1}{2}u_i)u_i - d_i x)dt \rightarrow \mathop{\max}\limits_{u_1,u_2 }$
Solution, first write down the Hamiltonion function:
$H(x,\psi,u)=\sum_{i=1}^{2}[ (b_i −\frac{1}{2}u_i)u_i − dx+ \psi(u_i)] $
where $\psi$ is the adjoint variable and $d=d_1+d_2$
Taking the first derivative with respect to $u_i$, we get the expressions for the optimal controls:
$u_{i}^{*}=b_i+\psi,i=1,2.$
Since there have no terminal cost in this case, then $\psi(T)=0.$ Combined with $\dot \psi=-\frac{\partial H}{\partial x}=d$, denote $d=d_1+d_2, b=b_1+b_2,$ therefore,
$\psi(t)=d(t-T).$
Correspondingly, the optimal control is
$u^{*}(t)=b_i-d(T-t)$, there still have one additional condition,
$d_i \in [0,\frac{min\left \{ b_1,b_2 \right \}}{T}-d]$,
The initial condition $x(0)=x_0$, naturally
$x^{*}(t)= \frac{2d}{2}(t^{2}-t_{0}^{2})+(b-2Td)(t-t_0)+x_0.$
Case 2， Consider a similar pollution control model. There are 2 players join in the game, N = {1, 2}. Each player has an industrial production site. It is assumed that the production is proportional to the pollutions $v_i$. Thus, the strategy of a player is to choose the amount of pollutions emited to the atmosphere, $v_i ∈ [0; b_i]$. In this example the solution will be considered in the class of open-loop strategies $v_i(t)$.
The dynamics of the total amount of pollution x(t) is described by
$\dot x=v_1+v_2-\delta x, x(t_0)=x_0.$
The payoff of the i-th player is defined as
$\int_{0}^{T} (p_i −\frac{1}{2}v_i)v_i − g_ix)dt, i = 1, 2.$
We assume that the terminal cost is zero.
Each player is to maximize their total cost, then the optimization problem is as follows:
$\sum_{i=1}^{2}\int_{0}^{T}(p_i-\frac{1}{2}v_i)v_i - g_i x)dt \rightarrow \mathop{\max}\limits_{v_1,v_2 }$
Similar to above, we use Pontrygin's maximal principle to get
$v_i^{*}(t)=p_i-\frac{g}{\delta}+\frac{g}{\delta}e^{\delta(t-T)},i=1,2$
Here, we assume that for the environment's self cleaning capacity $\delta$ the following condition hold: $\delta\geq \frac{d}{\ b_i}$. Then the optimal cooperative trajectory is
$x^{*}(t)= Ce^{-\delta t} + \frac{nd}{2\delta^2}e^{\delta(t-T)}+\frac{p}{\delta}-\frac{2g}{\delta^2}$
where $C=e^{\delta t_0}(x_0-\frac{p}{\delta}+\frac{2g}{\delta^2}-\frac{g}{\delta^2}e^{\delta(t_0-T)}), p=p_1+p_2, g=g_1+g_2$.

Question: If I combine two different model;
Dynamic system
$\dot x=u_1+u_2, x(t_0)=x_0.$
$\dot y=v_1+v_2-\delta y, y(t_0)=y_0$.
Optimization problem as follows:
$\rightarrow \mathop{\max}\limits_{u_1,u_2,v_1,v_2 }\sum_{i=1}^{2}\int_{0}^{T}[(p_i-\frac{1}{2}v_i)v_i - g_i y+(b_i −\frac{1}{2}u_i)u_i − d_i x)]dt $
where $p_i,g_i, b_i, d_i$ are parameters
We assume the there doesn't have terminal cost, and how to solve this kind of problem, previous I guess it is similar to finite Horizon LQ Game, but detail was different, does anyone can be give me some hints and advice?
 A: $$
\def\eqdef{\stackrel{\text{def}}{=}}
$$
First note that your cost functional doesn't depend on the individual values of $\ g_i\ $ and $\ d_i\ $, but only on the sums $\ g\eqdef g_1+g_2\ $ and $\ d\eqdef d_1+d_2\ $.
By integrating the differential equations for $\ x\ $ and $\ y\ $,
\begin{align}
x(t)&=x_0+\int_{t_0}^t\big(u_1(s)+u_2(s)\big)\,ds\\
y(t)&=e^{\delta\big(t_0-t\big)}y_0+\int_{t_0}^t e^{\delta(s-t)}\big(v_1(s)+v_2(s)\big)\,ds\ ,
\end{align}
and substituting the results into the cost functional, the latter can be written as
\begin{align}
C(u,v)&=C_0+\int_{t_0}^T\sum_{i=1}^2\left(\xi_iv_i+\phi_iu_i-\frac{1}{2}\big(v_i^2+u_i^2\big)\right)\,dt\\
&=C_1-\frac{1}{2}\int_{t_0}^T\sum_{i=1}^2\big(\big(v_i-\xi_i\big)^2+\big(u_i-\phi_i\big)^2\big)\,dt\ ,
\end{align}
where
\begin{align}
C_0&\eqdef\frac{gy_0}{\delta}\left(1-e^{t_0-T}\right)+dx_0\big(T-t_0\big)\ ,\\
C_1&\eqdef C_0+\frac{1}{2}\int_{t_0}^T\sum_{i=1}^2\big(\xi_i^2+\phi_i^2\big)\,dt\ ,\\
\xi_i(t)&\eqdef p_i-\frac{g}{\delta}\left(1-e^{\delta(t-T)}\right)\ \text{, and}\\
\phi_i(t)&\eqdef b_i-d(T-t)\ .
\end{align}
It's easy to show that the maximum of the cost function is achieved when
\begin{align}
v_i(t)&=\cases{0&if $\ \xi_i(t)<0$\\
               \xi_i(t)&if $\ 0\le\xi_i(t)\le p_i$\\
    p_i&if $\ p_i<\xi_i(t)$}\\
    \\
 u_i(t)&=\cases{0&if $\ \phi_i(t)<0$\\
               \phi_i(t)&if $\ 0\le\phi_i(t)\le b_i$\\
    p_i&if $\ b_i<\phi_i(t)\ $.}
\end{align}
If $\ \frac{g}{\delta}>0\ $ and $\ d>0\ $ then $\ \xi_i\ $, $\ \phi_i\ $ are all strictly increasing functions with $\ \xi_i(T)=p_i\ $, $ \ \phi_i(T)=b_i\ $, and $\ \phi_i\big(\tau_{u_i}\big)=0\ $, where
\begin{align}
\tau_{u_i}\eqdef T-\frac{b_i}{d}\ .
\end{align}
If $\ \frac{g}{\delta}\le p_i\ $, then $\ 0\le\xi_i(t)\le p_i\ $  for all $\ t\le T\ $ and the optimal values of $\ v_i, u_i\ $ are given by
\begin{align}
v_i(t)&=p_i-\frac{g}{\delta}\left(1-e^{\delta(t-T)}\right)
    \\
 u_i(t)&=\cases{0&if $\ t_0\le t<\tau_{u_i}$\\
    b_i-d(T-t)&if $\ \max\big(t_0,\tau_{u_i}\big)\le t\ $.}
\end{align}
If $\ \frac{g}{\delta}> p_i\ $, then $\ 
\xi_i\big(\tau_{v_i}\big)=0\ $, where
\begin{align}
\tau_{v_i}&\eqdef T-\frac{1}{\delta}\ln\left(1-\frac{\delta p_i}{g}\right)\  ,
\end{align}
and when $\ t_0<\tau_{v_i}\ $ the optimal value of $\ v_i(t)\ $ will be zero for $\ t_0\le t\le\tau_{v_i}\ $.
Reply to OP'S query in comments below (too long for a comment)
While I'm sure you could use both Pontryagin's maximum principle and the Hamilton-Jacobi-Bellman equation to derive the solution of your problem, I didn't use either of them. Nor did I use a computer program.  Your problem is a special case of one of the form
\begin{align}
&\max_\omega\hspace{3em}c_0-\|\ell-\omega\|^2\\
&\text{subject to}\hspace{1em} \omega\in\Omega\ ,
\end{align}
where $\ \Omega\ $ is a closed, strictly convex subset of the real Hilbert space $\ \mathscr{H}\eqdef\mathscr{L}^2\big([0,T],\mathbb{R}^n\big)\ $ of functions $\ f:[0,T]\rightarrow\mathbb{R}^n\ $ with $\ \int_0^T\|f(t)\|^2\,dt<\infty\ $, and $\ \ell\ $ some given element of $\ \mathscr{H}\  $.  Under these conditions there's a unique $\ w^*\ $ of $\ \Omega\ $ which is closest to $\ \ell\ $, and this will also be the unique member of $\ \Omega\ $ that maximises the objective function.  In your case, $\ n=4\ $, and the set $\ \Omega\ $ has the form
\begin{align}\big\{\omega\in\mathscr{H}\,\big|\,0\le\omega_i(t)\le\beta_i\ \text{ for }\ 0\le t\le T, i=1,2,\dots,n\big\}\ .
\end{align}
For this problem it's easy to show that $\ \omega^*\ $ is given by
$$
\omega_i^*(t)=\cases{0&if $\ \ell_i(t)<0$\\
       \ell_i(t)&if $\ 0\le\ell_i(t)\le\beta_i$\\
   \beta_i&if $\ \beta_i<\ell_i(t)$ . }
$$
