# Non-Markovian Reinforcement Learning using Fractional Dynamics

Gaurav Gupta  Ming Hsieh Department of Electrical and Computer Engineering,
University of Southern California, Los Angeles, CA, USA  Chenzhong Yin  Ming
Hsieh Department of Electrical and Computer Engineering,
University of Southern California, Los Angeles, CA, USA  Jyotirmoy V. Deshmukh
Department of Computer Science, University of Southern California, Los
Angeles, CA, USA  Paul Bogdan  Ming Hsieh Department of Electrical and
Computer Engineering,
University of Southern California, Los Angeles, CA, USA

###### Abstract

Reinforcement learning (RL) is a technique to learn the control policy for an
agent that interacts with a stochastic environment. In any given state, the
agent takes some action, and the environment determines the probability
distribution over the next state as well as gives the agent some reward. Most
RL algorithms typically assume that the environment satisfies Markov
assumptions (i.e. the probability distribution over the next state depends
only on the current state). In this paper, we propose a model-based RL
technique for a system that has non-Markovian dynamics. Such environments are
common in many real-world applications such as in human physiology, biological
systems, material science, and population dynamics. Model-based RL (MBRL)
techniques typically try to simultaneously learn a model of the environment
from the data, as well as try to identify an optimal policy for the learned
model. We propose a technique where the non-Markovianity of the system is
modeled through a fractional dynamical system. We show that we can quantify
the difference in the performance of an MBRL algorithm that uses bounded
horizon model predictive control from the optimal policy. Finally, we
demonstrate our proposed framework on a pharmacokinetic model of human blood
glucose dynamics and show that our fractional models can capture distant
correlations on real-world datasets.

##  1 Introduction

Reinforcement learning (RL) [1] is a technique to synthesize control policies
for autonomous agents that interact with a stochastic environment. The RL
paradigm now contains a number of different kinds of algorithms, and has been
successfully used across a diverse set of applications including autonomous
vehicles, resource management in computer clusters [2], traffic light control
[3], web system configuration [4], and personalized recommendations [5]. In
RL, we assume that in each state, the agent performs some action and the
environment picks a probability distribution over the next state and assigns a
reward (or negative cost). The reward is typically defined by the user with
the help of a state-based (or state-action-based) reward function. The
expected payoff that the agent may receive in any state can be defined in a
number of different ways; in this paper, we assume that the payoff is an
discounted sum of the local rewards (with some discount factor
$\gamma\in[0,1]$) over some time horizon $H$. The purpose of RL is to find the
stochastic policy (i.e. a distribution over actions conditioned on the current
state), that optimizes the expected payoff for the agent. Most RL algorithms
assume that the environment satisfies Markov assumptions, i.e. the probability
distribution over the next state is dependent only on the current state (and
not the history). In contrast, here, we investigate an RL procedure for a non-
Markovian environment.

Broadly speaking, there are two classes of RL algorithms [6]: model-based and
model-free algorithms. Most classical RL algorithms are model-based; they
assume that the environment is explicitly specified as a Markov Decision
Process (MDP), and use dynamic programming to compute the expected payoff for
each state of the MDP (called its value), as well as the optimal policy [7,
8]. Classical RL algorithms have strong convergence guarantees stemming from
the fact that the value of a state can be recursively expressed in terms of
the value of the next state (called the Bellman equation), which allows us to
define an operator to update the value (or the policy) for a given state
across iterations. This operator (also known as the Bellman operator) can be
shown to be a contraction mapping [1]. However, obtaining exact symbolic
descriptions of models is often infeasible. This led to the development of
model-free reinforcement learning (MFRL) approaches that rely on sampling many
model behaviors through simulations and eschew building a model altogether.
MFRL algorithms can converge to an optimal policy under the right set of
assumptions; however, can suffer from high sample complexity (i.e. the number
of simulations required to learn an optimal policy). This has led to
investigation of a new class of model-based RL (MBRL) algorithms where the
purpose is to simultaneously learn the system model as well as the optimal
policy [2]. Such algorithms are called on policy, as the policy learned during
any iteration is used for improving the learned model as well as optimizing
the policy further. Most MBRL approaches use function approximators or
Bayesian models to efficiently learn from scarce sample sets of system
trajectories. MBRL approaches tend to have lower sample complexity than MFRL
as the learned model can accelerate the convergence by focusing on actions
that are likely to be close to the optimal action. However, MBRL approaches
can suffer severely from modeling errors [9], and may converge to less optimal
solutions.

Figure 1: Non-Markovian Model Based Reinforcement Learning setup. The model
based predictions are used to select actions, and then iteratively update the
model dynamics.

In both MFRL and MBRL algorithms, a fundamental assumption is that the
environment satisfies Markovian properties, partly to avoid the complexity of
dealing with the historical dependence in transitions. To overcome this
challenge, we propose a non-Markovian MBRL framework that captures non-
Markovian characteristics through a fractional dynamical systems formulation.
Fractional dynamical systems can model non-Markovian processes characterized
by a single fractal exponent and commonly arise in mathematical models of
human physiological processes [10, 11], biological systems, condensed matter
and material sciences, and population dynamics [12, 13, 14, 15]. Such systems
can effectively model spatio-temporal properties of physiological signals such
as blood oxygenation level dependent (BOLD), electromyogram (EMG),
electrocardiogram (ECG), etc. [12, 16, 17]. The advantage of using fractional
dynamical models is that they can accurately represent long-range (historical)
correlations (memory) through a minimum number of parameters (e.g., using a
single fractal exponent to encode a long-range historical dependence rather
than memorizing the trajectory itself or modeling it through a large set of
autoregressive parameters). Though fractional models can be used to perform
predictive control [18], problems such as learning these models effectively or
obtaining optimal policies for such models in an RL setting have not been
explored.

In this paper, we develop a novel non-Markovian MBRL technique in which our
algorithm alternates between incrementally learning the fractional exponent
from data and learning the optimal policy on the updated model. We show that
the optimal action in a given state can be efficiently computed by solving a
quadratic program over a bounded horizon rollout from the state. The overview
of our model-based reinforcement learning algorithm is shown in Fig. 1. In
this algorithm, we use on-policy simulations to gather additional RL data that
is then used to update the model. Our model learning algorithm is based on
minimizing the distance between the data’s state-action distribution and the
next state distribution induced by the controller. The fractional dynamic
model is then retrained using the cumulative dataset. The MBRL procedure is
run for a finite number of user-specified iterations.

The rest of this paper is constructed as follows. We present our problem
statement in Section 2. Section 3 contains our proposed non-Markovian MBRL
algorithm. We demonstrate our experimental results in Section 4. In the end,
we conclude this paper with discussion and conclusion in Section 5.

##  2 Problem Formulation

The reinforcement learning deals with the design of the controller (or policy)
which minimizes the expected total cost. In the setting of a memoryless
assumption, the Markov Decision Process (MDP) [19] is used to model the system
dynamics such that the future state depends only on the current state and
action. For a state ${\bf s}_{t}\in\mathbb{R}^{n}$ and action ${\bf
a}_{t}\in\mathbb{R}^{p}$, the future state evolve as ${\bf s}_{t+1}\sim P({\bf
s}_{t+1}|{\bf s}_{t},{\bf a}_{t})$, and a cost function $r_{t}=c({\bf
s}_{t},{\bf a}_{t})$. However, the Markov assumption does not work well with
the long-range memory processes [20]. In this work, we take the non-Markovian
setting, or History Dependent Process (HDP), and hence, the future state
depends not only on the current action but also the history of states. The
history at time $t$ is the set $\mathcal{H}_{t}=\\{(s_{k})_{k\leq t}\\}$, and
for a trajectory $h\in\mathcal{H}_{t}$, we have $P({\bf s}_{t+1}|h,{\bf
a}_{t})$, or alternatively, $P_{h}({\bf s}_{t+1}|s_{t},{\bf a}_{t})$, where
the terminal state of the trajectory $h$ is written as $h(t)=s_{t}$. We
consider a model-based approach for reinforcement learning in a finite-horizon
setting. A non-Markovian policy $\pi(.|h)$ provides a distribution over
actions given the history of states until time $t$ as $h\in\mathcal{H_{t}}$.
For a given policy, the value function is defined as
$V_{h}^{\pi}=\mathbb{E}_{\pi(.|h)}\sum\nolimits_{t=0}^{T-1}c(s_{t},a_{t})$,
where the expectation is taken over state trajectories using policy $\pi$ and
the HDP, and $T$ is the horizon under consideration. We formally define the
non-Markovian MBRL problem in the Section 2.2.

###  2.1 Fractional Dynamical Model

A linear discrete time fractional-order dynamical model is described as
follows:

$\Delta^{\alpha}{\bf s}[k+1]={\bf A}{\bf s}[k]+{\bf B}{\bf a}[k],$ (1)

where ${\bf s}\in\mathbb{R}^{n}$ is the state, ${\bf a}\in\mathbb{R}^{p}$ is
the input action. The difference between a classic linear time-invariant (or
Markovian) and the above model is the inclusion of fractional-order derivative
whose expansion and discretization for any $i$th state $(1\leq i\leq n)$ can
be written as

$\Delta^{\alpha_{i}}s_{i}[k]=\sum\limits_{j=0}^{k}\psi(\alpha_{i},j)s_{i}[k-j],$
(2)

where $\alpha_{i}$ is the fractional order corresponding to the $i$th state
dimension and
$\psi(\alpha_{i},j)=\frac{\Gamma(j-\alpha_{i})}{\Gamma(-\alpha_{i})\Gamma(j+1)}$
with $\Gamma(.)$ denoting the gamma function. The system dynamics can also be
written in the probabilistic manner as follows:

$\displaystyle P_{{\bm{\theta}}}({\bf s}[k+1]|{\bf s}[0],\ldots,{\bf
s}[k],{\bf a}[k])=\mathcal{N}({\bm{\mu}}_{{\bm{\theta}}},{\bm{\Sigma}}),$
$\displaystyle{\bm{\mu}}_{{\bm{\theta}}}=\begin{bmatrix}[l]\sum\nolimits_{j=1}^{k}\psi(\alpha_{i},j)s_{0}[k-j]+{\bf
a}_{0}^{T}{\bf s}[k]+{\bf b}_{0}^{T}{\bf a}[k]+\mu_{0}\\\
\sum\nolimits_{j=1}^{k}\psi(\alpha_{i},j)s_{1}[k-j]+{\bf a}_{1}^{T}{\bf
s}[k]+{\bf b}_{1}^{T}{\bf a}[k]+\mu_{1}\\\ \vdots\\\
\sum\nolimits_{j=1}^{k}\psi(\alpha_{i},j)s_{n-1}[k-j]+{\bf a}_{n-1}^{T}{\bf
s}[k]+{\bf b}_{n-1}^{T}{\bf a}[k]+\mu_{n-1}\end{bmatrix},$ (3)

where ${\bm{\theta}}=\\{\alpha,{\bf A},{\bf B},\mu,{\bm{\Sigma}}\\}$, and
${\bf A}=[{\bf a}_{0},\ldots,{\bf a}_{n-1}]$, ${\bf B}=[{\bf
b}_{0},\ldots,{\bf b}_{n-1}]$. The fractional differencing operator in (3)
introduce the non-Markovianity by having long-range filtering operation on the
state vectors.

###  2.2 Non-Markovian Model Based Reinforcement Learning

The actions in MBRL are preferred on the basis of predictions made by the
undertaken model of the system dynamics. For many real-world systems, example
blood glucose [18, 21], ECG activities [11], the assumption of Markovian
dynamics does not hold and hence the predictions are not accurate, leading to
less rewarding actions selected for the system. As we note in the previous
section 2.1 that non-Markovian dynamics can be effectively and compactly
modeled as fractional dynamical system, we aim to use this system model for
making predictions. The non-Markovian MBRL problem is formally defined as
follows.

Problem Statement: Given non-Markovian state transitions, and actions dataset
in the time horizon $k\in[0,T-1]$ as $\mathcal{D}=\big{\\{}({\bf
s}[0],\ldots,{\bf s}[k],{\bf a}[k]),{\bf s}[k+1]\big{\\}}$. Let
$P_{{\bm{\theta}}}({\bf s}[k+1]|{\bf s}[0],\ldots,{\bf s}[k],{\bf a}[k])$ be
the non-Markovian system dynamics parameterized by the model parameters
${\bm{\theta}}$. Estimate the optimal policy which minimizes the expected
future discounted cost

$\displaystyle\pi^{\ast}=\operatorname*{arg\,min}\limits_{\pi}\,\mathbb{E}\sum\limits_{k=0}^{T-1}\gamma^{k}c({\bf
s}[k],{\bf a}[k]),$ (4)

where $\gamma$ is the discount factor satisfying $\gamma\in[0,1]$, and T is
the horizon under consideration.

##  3 Non-Markovian Reinforcement Learning

The MBRL comprises of two key steps, namely (_i_) the estimation of the model
dynamics from the given data $\mathcal{D}$, and (_ii_) the design of a policy
for optimal action selection which minimizes the total expected cost using
estimated dynamics. We discuss the solution to the non-Markovian MBRL as
follows.

###  3.1 Non-Markovian Model Predictive Control

The Model Predictive Control (MPC) aims at estimating the closed-loop policy
by optimizing the future discounted cost under a limited-horizon $H$ using
some approximation of the environment dynamics and the cost. In this work, we
are concerned with HDP using non-Markovian state dynamics. In MPC, the policy
could be a deterministic action, or a distribution over actions, and we sample
the action at each time-step in the latter. The MPC problem to estimate the
policy at time-step $k$ for a given $h\in\mathcal{H}_{k}$ can be formally
defined as

$\displaystyle\begin{aligned}
&\min\limits_{\pi(.|h)}&&\sum\nolimits_{l=k}^{k+H-1}\gamma^{l-k}\hat{c}(s[l],a[l])\\\
&\text{subject to}&\\\ &&&s[l+1]=f(h,a[l],e[l]),\forall l\geq k\\\
\end{aligned}$ (5)

The approximation of the environment dynamics $f$ could be non-linear in
general, and $e[l]$ is the system perturbation noise following some
distribution $e\sim g_{e}$. The presence of $e$ provides randomness in the
action sampling through policy, and the sampled action at each step is $a[k]$.
The performance of the non-Markovian MPC based policy is bounded within the
optimal policy using the following result.

######  Theorem 1.

Given an approximate HDP with
$||\hat{P}_{h^{\prime}}(s^{\prime}|s,a)-P_{h}(s^{\prime}|s,a)||_{1}\leq\mathcal{O}(t^{q})$,
$\forall h,h^{\prime}\in\mathcal{H}_{t}$ with $h(t)=h^{\prime}(t)=s$, and
$||c(s,a)-\hat{c}(s,a)||_{\infty}\leq\varepsilon$. The performance of the non-
Markovian MPC based policy $\hat{\pi}$ is related to the optimal policy
$\pi^{*}$ as

$\displaystyle||V_{h_{0}}^{\hat{\pi}}-V_{h_{0}}^{\pi^{*}}||_{\infty}$
$\displaystyle\leq
2\frac{1-\gamma^{H}}{1-\gamma}\left(\frac{c_{max}-c_{min}}{2}\right)H\mathcal{O}(T^{q})+2\varepsilon\frac{1-\gamma^{H}}{1-\gamma}\frac{1-\gamma^{T}}{1-\gamma},$
(6)

where, $h_{0}\in\mathcal{H}_{0}$ is the initial history given to the system.

The assumption of model approximation is critical here, and the error
increases if the exponent $q$ increases. For the MDP setting, the
approximation is taken as $q=0$. However, for a HDP with the history of length
$t$, we scale the approximation gap with $t$. The MPC horizon also plays a
role in the error bound, and the error increases for larger $H$.

The non-Markovian MPC could be computationally prohibitive (expensive) in the
general setting. Consequently, we now discuss the fractional dynamical MPC
approach which is non-Markovian but computationally tractable.

###  3.2 Fractional Model Predictive Control

The linear discrete fractional dynamical model as discussed in (1) is used as
an approximation to the non-Markovian environment dynamics. Formally, for our
purpose, the fractional MPC problem using (5) is defined as

$\displaystyle\begin{aligned} &&\min\limits_{{\bf
a}[k]}\sum\nolimits_{l=k}^{k+H-1}&\gamma^{l-k}\hat{c}({\bf s}[l],{\bf
a}[l])\\\ &\text{s.t.}&\\\ &&\Delta^{\alpha}\bar{\bf s}[l+1]&={\bf A}\bar{\bf
s}[l]+{\bf B}{\bf a}[l]+e[l],\\\ &&\bar{\bf s}[k^{\prime}]&={\bf s}[k],\forall
k^{\prime}\leq k,\\\ &&{\bf s}_{min}&\leq\bar{\bf s}[l]\leq{\bf
s}_{max},\forall l,\end{aligned}$ (7)

where ${\bf s}_{min},{\bf s}_{max}$ are feasibility bounds on the problem
according to the application, and the model noise
$e\sim\mathcal{N}(0,\Sigma)$. Note that (7) provides a policy using fractional
MPC. The action $a[k]$ is sampled from this policy by first sampling
$e\sim\mathcal{N}(0,\Sigma)$, and then solving (7). The non-Markovian
fractional dynamics would introduce the computation complexities in optimally
solving the problem in (7). However, since the constraints in (7) are linear,
for cost approximations $\hat{c}$ that are quadratic, a quadratic programming
(QP) solution can be developed to solve the fractional MPC efficiently. We
refer the reader to Appendix B for the QP version of the fractional MPC.
Further, a convex formulation of the costs $\hat{c}$ also enables efficient
solution of the fractional MPC using convex programming solvers, for example,
CPLEX and Gurobi [22, 23].

Next, we discuss the methodologies required to make an approximation of the
non-Markovian environment using fractional dynamics.

###  3.3 Model Estimation

The fractional dynamical model as described in the Section 2.1 is estimated
using the approach proposed in [12] by replacing the unknown inputs with known
actions at any time-step. For the sake of completeness, we present estimation
algorithm as Algorithm 1. We note that in [12] the input data is obtained only
once, and hence in this work appropriate modification in Algorithm 1 is
performed to work with recursively updated dataset as we see in Section 3.4.

1Input: $\mathcal{D}=\big{\\{}({\bf s}[0],\ldots,{\bf s}[k],{\bf a}[k]),{\bf
s}[k+1]\big{\\}}$ in the time-horizon $k\in[0,T-1]$

Output: ${\bm{\theta}}=\\{\alpha,{\bf A},{\bf B},\mu,{\bf\Sigma}\\}$

1: Estimate $\alpha$ using wavelets fitting for each state dimension

2: for $i=1,2,\ldots,n$ do

3: Compute $z_{i}[k]=\Delta^{\alpha_{i}}s[k+1]$ using $\alpha_{i}$
$\triangleright$ Eq.(2)

4: Aggregate $z_{i}[k],s[k],a[k]$ as $Z_{i},S,U$

5: $[a_{i}^{T},b_{i}^{T},\mu]=\arg\min\limits_{a,b,\mu}||Z_{i}-Sa-
Ub-\mu||_{2}^{2}$ with ${\bm{\Sigma}}$ as squared error

6: end for

Algorithm 1 Fractional_Dynamics_Estimation

The Markovian model assume memoryless property and hence lacks long-range
correlations for further accurate modeling. The existence of long-range
correlations can be estimated by computing the Hurst exponent $\bar{H}$. For
long-range correlations, the $\bar{H}$ lies in the range of $(0.5,1]$. The
fractional coefficient $\alpha$ in our model is related with $\bar{H}$ as
$\alpha=\bar{H}-0.5$. The Hurst exponent can be estimated from the slope of
log-log variations of the variance of wavelets coefficients vs scale as noted
in [24]. In the experiments Section 4.2, we show log-log plot to observe the
presence of long-range correlations in the real-world data.

###  3.4 Model Based Reinforcement Learning

The non-Markovian MPC exploiting the fractional dynamical model formulation in
Section 3.2 utilizes a dataset of the form $\mathcal{D}=\big{\\{}({\bf
s}[0],\ldots,{\bf s}[k],{\bf a}[k]),{\bf s}[k+1]\big{\\}}$ in the time-horizon
$k\in[0,T-1]$. We note that the performance of such MPC can be further
improved by using reinforcement learning. The selected actions by the MPC
${\bf a}[k]$ can be used to gather new transitions ${\bf s}[k+1]|{\bf
s}[0],\ldots,{\bf s}[k],{\bf a}[k]$, or acquiring data using on-policy. The
aggregated data is now used to re-estimate the model dynamics, and then
perform MPC. Specifically, the MBRL proceeds as follows. Using the seed
dataset, a parameterized fractional model dynamics is estimated as
$P_{{\bm{\theta}}}({\bf s}[k+1]|{\bf s}[0],\ldots,{\bf s}[k],{\bf a}[k])$. The
model dynamics is used to minimize the discounted future cost as MPC in
equation (7). The selected action along with the history of states ${\bf
s}[0],\ldots,{\bf s}[k]$ is used to gather the next transition using on-policy
as ${\bf s}[k+1]|{\bf s}[0],\ldots,{\bf s}[k],{\bf a}[k]$. The seed dataset is
updated with the gathered on-policy data $\mathcal{D}_{RL}$ to get aggregated
dataset. The fractional dynamics are updated using the new dataset, and the
aforementioned steps are repeated for a given number of iterations. The above
steps are summarized as Algorithm 2. The Algorithm 2 utilizes Algorithm 1
iteratively for the fractional model estimation. We now proceed to Section 4
for numerical demonstration of the proposed schemes.

1Input: Seed dataset $\mathcal{D}_{s}=\big{\\{}({\bf s}[0],\ldots,{\bf
s}[k],{\bf a}[k]),{\bf s}[k+1]\big{\\}}$ in the time-horizon $k\in[0,T-1]$

2 Output: ${\bm{\theta}}$

Initialize: $\mathcal{D}_{RL}\leftarrow\phi$

1: for $iter=1,2,\ldots,iter\\_max$ do

2:
$\theta\leftarrow\text{Fractional}\\_\text{Dynamics}\\_\text{Estimation}(\mathcal{D}_{s}\cup\mathcal{D}_{RL})$

3: Set initial state $\bar{\bf s}[0]\leftarrow{\bf s}[0]$

4: for $k=0,1,\ldots,T-1$ do

5: Sample action $a[k]$ from the fractional MPC based policy using $\bar{{\bf
s}}[l],\forall l\leq k$ $\triangleright$ Eq.(7)

6: Get $\bar{\bf s}[k+1]$ by executing ${\bf a}[k]$

7: $\mathcal{D}_{RL}\leftarrow\mathcal{D}_{RL}\cup\big{\\{}(\bar{\bf
s}[0],\ldots,\bar{\bf s}[k],{\bf a}[k]),\bar{\bf s}[k+1]\big{\\}}$

8: end for

9: end for

Algorithm 2 Fractional_Reinforcement_Learning

##  4 Experiments

We show example of the fractional MBRL on the blood glucose (BG) control. The
motive of blood glucose control is to make the BG in the range of
$70-180mg/dL$. The BG control is crucial in the treatment of T1 diabetes
patients which have inability to produce the required insulin amounts. The low
levels of glucose in the blood plasma is termed as hypoglycemia, while the
high levels is termed as hyperglycemia. For the application of reinforcement
learning, the cost function is taken as risk associated with different levels
of BG in the system. In [25] a quantified version of risk is proposed as
function of BG levels which is written as follows.

$\displaystyle f(b)$ $\displaystyle=1.509\times(\log(b)^{1.084}-5.381),$
$\displaystyle R(b)$ $\displaystyle=10\times(f(b))^{2}.$ (8)

Next, the cost for the transition instance $s[k+1]|s[0],\ldots,s[k],a[k])$ is
written as

$\hat{c}(s[k],a[k])=R(s[k+1])-R(s[k]),$ (9)

where the state $s[k]\in\mathbb{R}$ represents the BG level at time instant
$k$, and $a[k]$ represents the insulin dose and $R(.)$ is from (8). In rest of
the section, we experiment with simulated and real-world dataset,
respectively.

###  4.1 UVa T1DM Simulator

$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $76.67\%$ | $23.32\%$
$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $84.52\%$ | $15.47\%$
$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $93.53\%$ | $6.46\%$

(a)

$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $74.13\%$ | $25.86\%$
$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $87.52\%$ | $12.47\%$
$<70$ | $70-180$ | $>180$
---|---|---
$mg/dL$ | $mg/dL$ | $mg/dL$
$0\%$ | $90.30\%$ | $9.69\%$

(b)

Figure 2: Blood Glucose (BG) level with time, by implementation of fractional
Reinforcement Learning Scheme as Controller, of two Adults in (a) and (b). For
each subject, the BG level trajectories are shown from left-to-right in the
increasing number of RL iterations with leftmost, middle, and rightmost are
outputs at $5$, $10$, and $15$ iterations. As RL iterations increase the MBRL
scheme learns better policy and the BG level stays more in the desired level
of $70-180mg/dL$. The percentage of time spend in different BG level zone is
shown in tables below each plot.

The UVa/Padova T1DM [26] is a FDA approved T1 Diabetes simulator which
supports multiple virtual subjects. An open-source implementation of the
simulator [27] is used in this work. We take similar simulation setup as in
[28]. Each subject is simulated for a total of $36$ hours starting from 6 a.m.
in the morning. The meal timings/quantity are fixed as $50g$ CHO at 9 a.m.,
$70g$ at 1 p.m, $90g$ at 5:30 p.m, and $25g$ at 8 p.m. On day 2, $50g$ at 9
a.m., and $70g$ at 1 p.m. The continuous glucose monitor (CGM) sensor measures
the BG at every 5 mins.

For applying Algorithm 2, we set the horizon length $H$ in MPC be $100$
samples, discount factor $\gamma=0.99$. The $s_{min},s_{max}$ in MPC problem
(7) are set as $70,180$ respectively. The maximum number of RL iterations
$iter\\_max$ are set as $30$. We show the BG output of the simulator using
Algorithm 2 as controller in Fig. 2. We observe that the fraction of time BG
stays in the desired zone $70-180mg/dL$ increase with increasing the learning
iterations in the Algorithm 2. The data gathered using on-policy helps the
model making better prediction, and with as few iterations as $15$ we have
more than $90\%$ of time BG stays in the desired levels.

###  4.2 Real-World Data

Testing the controllers on real-world systems is difficult because of the
health risks associated with the patients. We take the Diabetes dataset from
UCI repository [29] which records the BG level and insulin dosage for $70$
patients. While testing controller is not possible here, hence we present the
analysis regarding the modeling part. The long-range memory in the signals
exist if the associated fractional exponent lies in the range of $(0,0.5]$ as
noted in Section 3.3. In Fig. 3, we show the log-log plots of the variance of
wavelets coefficients at various scales, for two subjects. We observe that the
estimated value of $\alpha$ lies in $(0,0.5]$ which indicates presence of
long-range memory, and hence fractional models can be used to make better
predictions.

###  4.3 Discussion

Insulin dependent diabetes mellitus (IDDM) is a kind of chronic disease
characterized by abnormal BG level. Topically, high level of BG, which is
caused by either the pancreas does not compound enough insulin (a hormone that
signals cells to uptake glucose in the bloodstream) or the produced insulin
cannot be effectively used by the human body, can result in a disorder of
metabolic that give rise to irreversible damage (such as lesion of patients’
organs, retinopathy, nephropathy, peripheral neuropathy and blindness) [30,
31]. According to recent research, nowadays, IDDM is influencing 20-40 million
people around the world and this amount is increasing over time [32]. Related
work [33] presents that tight blood glucose control along with insulin
injections can help to control the disease, however intensive control can
result in the risk of low blood sugar. This symptom can increase the risk of
heart disease, or even sudden death.

Figure 3: The log-log plot of variance of wavelet coefficients vs scale of two
subjects in (a) and (b). The values of $\alpha$ lies in $(0,0.5]$ which
indicates long-range correlations.

To effectively and safely against with IDDM, in the work of [34], the authors
constructed a CGM to detect the insulin among in the individuals in real time.
This monitor can read the blood glucose of the patients for every 5 mins.
Combining with an insulin pump (a small device that automatically inject
insulin), the CGM can constructed a system called "artificial pancreas" (AP).
The AP system is designed to control the symptom of patients which can
dynamically predict the among of insulin the individual should be delivered.
For many years, researchers have worked on designing efficient
algorithms/models to correctly predict the required insulin of individuals in
AP systems. In this paper, we present an innovative MBRL algorithm which is
explored in the non-Markovian model to dynamically make the prediction of the
among of insulin patients demand with high-accuracy and high-efficiency.

##  5 Conclusion

There are many important learning control problems that are not naturally
formulated as Markov decision processes. For example, if the agent cannot
directly observe the environment state, then the use of a partially observable
Markov decision process (POMDP) [35] model is more appropriate. Even in
presence of full observability, the probability distribution over next states
may not depend only on the current state. A more general class can be termed
as History Dependent Process (HDP), which can be looked as infinite-state
POMDP [36]. Another non-MDP class for model-free is Q-value Uniform Decision
Process (QDP) [37]. The non-Markovianity in the rewards structure is explored
in [38, 39] which utilize model-free learning, and RL for POMDP is explored in
[40] which is also model-free. MBRL is used for various robotics application
[41] in the MDP setting. The deep probabilistic networks using MDP is used in
[6].

In this work, we constructed a non-Markovian Model Based Reinforcement
Learning (MBRL) algorithm consisted with fractional dynamics model and the
model predictive control. The current Reinforcement learning (RL) approaches
have two kinds of limitations: (i) model-free RL models can achieve a high
predict accuracy, but these approaches need a large number of data-points to
train the model; (ii) current models don’t make latent behavioral patterns
into considerations which can affect the prediction accuracy in MBRL. We show
that our non-Markovian MBRL model can validly avoid these limitations.
Firstly, in our algorithm, we gather additional on-policy data to alternate
between gathering the initial data, hence it needs less sample points than the
general model-free RL approaches. Secondly, fractional dynamical model is the
key element in our algorithm to improve/guarantee the prediction accuracy. The
experiments on the blood glucose (BG) control to dynamically predict the
desired insulin amount show that the proposed non-Markovian framework helps in
achieving desired levels of BG for longer times with consistency.

The richness of complex systems cannot be always modeled as Markovian
dynamics. Previous works have shown that the long-range memory property of
fractional differentiation operators can model biological signals
efficaciously and accurately. Thus, we have modeled the blood glucose as non-
Markovian fractional dynamical system and developed solutions using
reinforcement learning approach. Finally, while the application of non-
Markovian MBRL open venues for real-world implementation but proper care has
to be taken especially when we have to deal with the healthcare systems. The
future investigations would involve more personalized modeling capabilities
for such systems with utilization of the domain knowledge. Nonetheless, we
show that the use of long-range dependence in the biological models is worth
exploring and simple models yield benefits of compactness as well as better
accuracy of the predictions.

## References

  * [1] D. P. Bertsekas, _Reinforcement learning and optimal control_. Athena Scientific Belmont, MA, 2019.
  * [2] H. Mao, M. Alizadeh, I. Menache, and S. Kandula, “Resource management with deep reinforcement learning,” in _Proceedings of the 15th ACM Workshop on Hot Topics in Networks_ , 2016, pp. 50–56.
  * [3] I. Arel, C. Liu, T. Urbanik, and A. G. Kohls, “Reinforcement learning-based multi-agent system for network traffic signal control,” _IET Intelligent Transport Systems_ , vol. 4, no. 2, pp. 128–135, 2010.
  * [4] X. Bu, J. Rao, and C.-Z. Xu, “A reinforcement learning approach to online web systems auto-configuration,” in _2009 29th IEEE International Conference on Distributed Computing Systems_. IEEE, 2009, pp. 2–11.
  * [5] G. Zheng, F. Zhang, Z. Zheng, Y. Xiang, N. J. Yuan, X. Xie, and Z. Li, “Drn: A deep reinforcement learning framework for news recommendation,” in _Proceedings of the 2018 World Wide Web Conference_ , 2018, pp. 167–176.
  * [6] K. Chua, R. Calandra, R. McAllister, and S. Levine, “Deep reinforcement learning in a handful of trials using probabilistic dynamics models,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 4754–4765.
  * [7] R. S. Sutton, D. Precup, and S. Singh, “Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning,” _Artificial intelligence_ , vol. 112, no. 1-2, pp. 181–211, 1999.
  * [8] A. L. Strehl, L. Li, and M. L. Littman, “Reinforcement learning in finite mdps: Pac analysis.” _Journal of Machine Learning Research_ , vol. 10, no. 11, 2009.
  * [9] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 5026–5033.
  * [10] B. West, “Fractal physiology and the fractional calculus: A perspective,” _Frontiers in Physiology_ , vol. 1, p. 12, 2010. [Online]. Available: https://www.frontiersin.org/article/10.3389/fphys.2010.00012
  * [11] Y. Xue, S. Pequito, J. R. Coelho, P. Bogdan, and G. J. Pappas, “Minimum number of sensors to ensure observability of physiological systems: a case study,” in _Allerton_ , 2016.
  * [12] G. Gupta, S. Pequito, and P. Bogdan, “Dealing with unknown unknowns: Identification and selection of minimal sensing for fractional dynamics with unknown inputs,” _in American Control Conference_ , 2018, arXiv:1803.04866.
  * [13] C. Yin, G. Gupta, and P. Bogdan, “Discovering laws from observations: A data-driven approach,” in _International Conference on Dynamic Data Driven Application Systems_. Springer, 2020, pp. 302–310.
  * [14] G. Gupta, S. Pequito, and P. Bogdan, “Learning latent fractional dynamics with unknown unknowns,” in _2019 American Control Conference (ACC)_ , 2019, pp. 217–222.
  * [15] ——, “Re-thinking eeg-based non-invasive brain interfaces: modeling and analysis,” in _2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS)_. IEEE, 2018, pp. 275–286.
  * [16] D. Baleanu, J. A. T. Machado, and A. C. Luo, _Fractional dynamics and control_. Springer Science & Business Media, 2011.
  * [17] R. L. Magin, _Fractional calculus in bioengineering_. Begell House Redding, 2006, vol. 2, no. 6.
  * [18] M. Ghorbani and P. Bogdan, “Reducing risk of closed loop control of blood glucose in artificial pancreas using fractional calculus,” in _2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , 2014, pp. 4839–4842.
  * [19] R. BELLMAN, “A markovian decision process,” _Journal of Mathematics and Mechanics_ , vol. 6, no. 5, pp. 679–684, 1957.
  * [20] S. Micciche, “Modeling long-range memory with stationary markovian processes,” _Physical Review E_ , vol. 79, no. 3, p. 031116, 2009.
  * [21] M. Otoom, H. Alshraideh, H. M. Almasaeid, D. López-de Ipiña, and J. Bravo, “A real-time insulin injection system,” in _International Workshop on Ambient Assisted Living_. Springer, 2013, pp. 120–127.
  * [22] J. Kronqvist, D. E. Bernal, A. Lundell, and I. E. Grossmann, “A review and comparison of solvers for convex minlp,” _Optimization and Engineering_ , vol. 20, no. 2, pp. 397–455, 2019.
  * [23] F. Hutter, H. H. Hoos, and K. Leyton-Brown, “Automated configuration of mixed integer programming solvers,” in _International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming_. Springer, 2010, pp. 186–202.
  * [24] P. Flandrin, “Wavelet analysis and synthesis of fractional brownian motion,” _IEEE Transactions on Information Theory_ , vol. 38, no. 2, pp. 910–917, March 1992.
  * [25] W. Clarke and B. Kovatchev, “Statistical tools to analyze continuous glucose monitor data,” _Diabetes Technol. Ther._ , vol. 11 Suppl 1, pp. 45–54, Jun 2009.
  * [26] B. P. Kovatchev, M. Breton, C. D. Man, and C. Cobelli, “In silico preclinical trials: a proof of concept in closed-loop control of type 1 diabetes,” _J Diabetes Sci Technol_ , vol. 3, no. 1, pp. 44–55, Jan 2009\.
  * [27] J. Xie, “Simglucose v0.2.1,” 2018. [Online]. Available: https://github.com/jxx123/simglucose
  * [28] Q. Wang, J. Xie, P. Molenaar, and J. Ulbrecht, “Model predictive control for type 1 diabetes based on personalized linear time-varying subject model consisting of both insulin and meal inputs: In silico evaluation,” in _2015 American Control Conference (ACC)_ , 2015, pp. 5782–5787.
  * [29] D. Dua and C. Graff, “UCI machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml
  * [30] M. Tejedor, A. Z. Woldaregay, and F. Godtliebsen, “Reinforcement learning application in diabetes blood glucose control: A systematic review,” _Artificial Intelligence in Medicine_ , p. 101836, 2020.
  * [31] M. Derouich and A. Boutayeb, “The effect of physical exercise on the dynamics of glucose and insulin,” _Journal of biomechanics_ , vol. 35, no. 7, pp. 911–917, 2002.
  * [32] W.-P. You and M. Henneberg, “Type 1 diabetes prevalence increasing globally and regionally: the role of natural selection and life expectancy at birth,” _BMJ open diabetes research and care_ , vol. 4, no. 1, 2016.
  * [33] D. Control, C. T. R. Group _et al._ , “Resource utilization and costs of care in the diabetes control and complications trial,” _Diabetes Care_ , vol. 18, no. 11, pp. 1468–1478, 1995.
  * [34] R. D. Coffen and L. M. Dahlquist, “Magnitude of type 1 diabetes self-management in youth health care needs diabetes educators,” _The Diabetes Educator_ , vol. 35, no. 2, pp. 302–308, 2009.
  * [35] R. D. Smallwood and E. J. Sondik, “The optimal control of partially observable markov processes over a finite horizon,” _Operations Research_ , vol. 21, no. 5, pp. 1071–1088, 1973.
  * [36] J. Leike, “Nonparametric general reinforcement learning,” 2016.
  * [37] S. J. Majeed and M. Hutter, “On q-learning convergence for non-markov decision processes,” in _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, IJCAI-18_. International Joint Conferences on Artificial Intelligence Organization, 7 2018, pp. 2546–2552.
  * [38] M. Gaon and R. I. Brafman, “Reinforcement learning with non-markovian rewards,” 2019.
  * [39] M. Agarwal and V. Aggarwal, “Reinforcement learning for joint optimization of multiple rewards,” 2021.
  * [40] J. Perez and T. Silander, “Non-markovian control with gated end-to-end memory policy networks,” 2017.
  * [41] A. Nagabandi, G. Kahn, R. S. Fearing, and S. Levine, “Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning,” _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2018. [Online]. Available: http://dx.doi.org/10.1109/ICRA.2018.8463189
  * [42] M. Bhardwaj, S. Choudhury, and B. Boots, “Blending mpc & value function approximation for efficient reinforcement learning,” 2020.
  * [43] S. Ross and J. A. Bagnell, “Agnostic system identification for model-based reinforcement learning,” 2012.
  * [44] S. Kakade, M. Kearns, and J. Langford, “Exploration in metric state spaces,” in _Proceedings of the Twentieth International Conference on International Conference on Machine Learning_ , ser. ICML’03, 2003, p. 306–312.

##  Appendix A Proof of Theorem 1

For the notational purpose of the proof, we define the non-Markovian
environment as $M$, and approximation to the environment as $\hat{M}$. The
non-Markovian MPC based policy is $\hat{\pi}$, and the optimal policy is
$\pi^{*}$. The approximation quality of the HDP dynamics are
$||\hat{P}_{h^{\prime}}(s^{\prime}|s,a)-P_{h}(s^{\prime}|s,a)||_{1}\leq\mathcal{O}(t^{q})$,
$\forall h,h^{\prime}\in\mathcal{H}_{t}$ with $h(t)=h^{\prime}(t)=s$, and the
approximation for cost $||c(s,a)-\hat{c}(s,a)||_{\infty}\leq\varepsilon$. We
also assume that the range of cost function is $[c_{min},c_{max}]$. The
initial history information provided is taken as $h_{0}\in\mathcal{H}_{0}$,
for example, initial state $s_{0}$. We define $\hat{h}(h^{*})$ as the
trajectories taken by policy $\hat{\pi}$ ($\pi^{*}$), respectively. We define
the value function at step-$t$ ($0\leq t\leq T-1$), with history h using the
model $M$, and policy $\pi$ as

$V_{t,M}^{\pi}(h)=\mathbb{E}_{h,\pi,M}\sum\limits_{k=t}^{T-1}\gamma^{k}(s_{k},a_{k}),$
(10)

where $h\in\mathcal{H}_{t}$ and $h(t)=s_{t}$. We first present the simulation
lemma for non-Markovian HDP as follows.

######  Lemma 1.

Given an approximation $\hat{M}$ of the environment $M$ as,
$||\hat{P}_{h^{\prime}}(s^{\prime}|s,a)-P_{h}(s^{\prime}|s,a)||_{1}\leq\mathcal{O}(t^{q})$,
$\forall h,h^{\prime}\in\mathcal{H}_{t}$ with $h(t)=h^{\prime}(t)=s$, and the
approximation for cost $||c(s,a)-\hat{c}(s,a)||_{\infty}\leq\varepsilon$, and
any policy $\pi$ and history $h_{t}\in\mathcal{H}_{t}$,

$\displaystyle||\mathbb{E}_{h,\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}(s_{k},a_{k})-\mathbb{E}_{h,\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}(s_{k},a_{k})||_{\infty}$
$\displaystyle\qquad\leq\gamma^{t}H\left(\frac{c_{max}-c_{min}}{2}\right)\frac{1-\gamma^{H}}{1-\gamma}\mathcal{O}((t+H)^{q})+\varepsilon\gamma^{t}\frac{1-\gamma^{H}}{1-\gamma}.$
(11)

###### Proof.

The difference in the value functions can be written as

$\displaystyle\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\qquad=\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\qquad+\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k}).$

With $||c(s,a)-\hat{c}(s,a)||_{\infty}\leq\varepsilon$, we can write that

$\displaystyle||\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\hat{c}(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}c(s_{k},a_{k})||_{\infty}$
$\displaystyle\qquad\leq\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\nolimits_{k=t}^{t+H-1}\gamma^{k}\varepsilon=\varepsilon\frac{1-\gamma^{H}}{1-\gamma}.$

For the remaining terms, we can write that

$\displaystyle\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\quad=\left(\sum\limits_{s_{k}\sim
h_{t},\pi,\hat{M}}{P}(s_{t},\ldots,s_{t+H-1})-\sum\limits_{s_{k}\sim
h_{t},\pi,M}{P}(s_{t},\ldots,s_{t+H-1})\right)(\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\delta),$

where in the last equality, we can insert $\delta$ as

$\displaystyle\sum\limits_{s_{k}\sim
h_{t},\pi,\hat{M}}{P}(s_{t},\ldots,s_{t+H-1})=\sum\limits_{s_{k}\sim
h_{t},\pi,M}{P}(s_{t},\ldots,s_{t+H-1})=0.$

Next,

$\displaystyle||\mathbb{E}_{h_{t},\pi,\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{h_{t},\pi,M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})||_{\infty}$
$\displaystyle\quad\leq||(\sum\limits_{s_{k}\sim
h_{t},\pi,\hat{M}}{P}(s_{t},\ldots,s_{t+H-1})-\sum\limits_{s_{k}\sim
h_{t},\pi,M}{P}(s_{t},\ldots,s_{t+H-1}))||_{\infty}\,\,||\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\delta||_{\infty}.$

By choosing
$\delta=\sum\nolimits_{k=t}^{t+H-1}\gamma^{k}(\frac{c_{max}+c_{min}}{2})$, and
upper-bounding difference of transition dynamics as $\mathcal{O}((t+H)^{q})$,
we get the final expression. ∎

Now, we begin the proof of the Theorem 1 as follows.

$\displaystyle
V_{t,M}^{\hat{\pi}}(\hat{h}_{t})-V_{t,M}^{\pi^{\ast}}(h^{*}_{t})$
$\displaystyle=\mathbb{E}_{\hat{h}_{t},\hat{\pi},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{h^{*}_{t},\pi^{*},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\qquad+V_{t+H,M}^{\hat{\pi}}(\hat{h}_{t+H})-V_{t+H,M}^{\pi^{\ast}}(h^{*}_{t+H}).$

The first set of terms can be expanded as follows.

$\displaystyle\mathbb{E}_{\hat{h}_{t},\hat{\pi},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{h^{*}_{t},\pi^{*},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\qquad=\mathbb{E}_{\hat{h}_{t},\hat{\pi},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{\hat{h}_{t},\hat{\pi},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})$
$\displaystyle\qquad\qquad+\mathbb{E}_{\hat{h}_{t},\hat{\pi},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})-\mathbb{E}_{\hat{h}_{t},\pi^{*},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})$
$\displaystyle\qquad\qquad+\mathbb{E}_{\hat{h}_{t},{\pi}^{*},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})-\mathbb{E}_{\hat{h}_{t},\pi^{*},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})$
$\displaystyle\qquad\qquad+\mathbb{E}_{\hat{h}_{t},{\pi}^{*},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k})-\mathbb{E}_{{h}^{*}_{t},\pi^{*},M}\sum\limits_{k=t}^{t+H-1}\gamma^{k}c(s_{k},a_{k}).$

Since the $\hat{\pi}$ is greedy policy that optimize (5), therefore,

$\mathbb{E}_{\hat{h}_{t},\hat{\pi},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k})\leq\mathbb{E}_{\hat{h}_{t},\pi^{*},\hat{M}}\sum\limits_{k=t}^{t+H-1}\gamma^{k}\hat{c}(s_{k},a_{k}).$
(12)

Now, using lemma 1 and the approximation quality of costs, we get

$\displaystyle||V_{t,M}^{\hat{\pi}}(\hat{h}_{t})-V_{t,M}^{\pi^{\ast}}(h^{*}_{t})||_{\infty}\leq$
$\displaystyle\quad
2\gamma^{t}H\left(\frac{c_{max}-c_{min}}{2}\right)\frac{1-\gamma^{H}}{1-\gamma}\mathcal{O}((t+H)^{q})+2\gamma^{t}\varepsilon\frac{1-\gamma^{H}}{1-\gamma}+||V_{t+H,M}^{\hat{\pi}}(\hat{h}_{t+H})-V_{t+H,M}^{\pi^{\ast}}(h^{*}_{t+H})||_{\infty}.$

By adding the terms from $t=0$ till $T-1$, the terms that are $H$ apart
cancels out in a telescopic sum fashion. Finally, we can write that

$\displaystyle||V_{0,M}^{\hat{\pi}}({h}_{0})-V_{0,M}^{\pi^{\ast}}(h_{0})||_{\infty}$
$\displaystyle\leq\sum\limits_{t=0}^{T-1}2\gamma^{t}H\left(\frac{c_{max}-c_{min}}{2}\right)\frac{1-\gamma^{H}}{1-\gamma}\mathcal{O}((t+H)^{q})+2\gamma^{t}\varepsilon\frac{1-\gamma^{H}}{1-\gamma}$
$\displaystyle\leq
2\frac{1-\gamma^{H}}{1-\gamma}\left(\frac{c_{max}-c_{min}}{2}\right)H\mathcal{O}(T^{q})$
$\displaystyle\qquad+2\varepsilon\frac{1-\gamma^{H}}{1-\gamma}\frac{1-\gamma^{T}}{1-\gamma}.$

∎

##  Appendix B Fractional MPC

The fractional MPC in Section 3.2 has linear constraints. Define the
optimization variable
$x=[\bar{s}[k+H]^{T},\bar{s}[k+H-1]^{T},\ldots,\bar{s}[0]^{T},a[k+H-1]^{T},\ldots,a[k]^{T}]^{T}$.
The first set of constraints can then be written as $\Theta\,x=b$, where

$\displaystyle\Theta=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}I&D(\alpha,1)+A&D(\alpha,2)&\ldots&&&&D(\alpha,k+H)&B&0&\ldots&0\\\
0&I&D(\alpha,1)+A&D(\alpha,2)&&&&\ldots&0&B&\ldots&0\\\
\vdots&\vdots&\ddots&\ldots&I&D(\alpha,1)+A&\ldots&D(\alpha,k+1)&0&0&\ldots&B\\\
\end{bmatrix},$

using,

$D(\alpha,j)=diag(\psi(\alpha_{1},j),\ldots,\psi(\alpha_{n},j)).$ (13)

and $b=[e[k+H-1]^{T},\ldots,e[k]^{T},0^{T},\ldots,0^{T}]^{T}$. The history
equality constraint can be set as $\Phi\,x=d$ with $\Phi=[0,I,0]$ of
appropriate size, and
$d=[0^{T},\ldots,s[k]^{T},\ldots,s[0]^{T},0^{T},\ldots,0^{T}]^{T}$. Using
these two equality constraints, and boundary limits for $\bar{s}[k]$, we get a
quadratic programming with approximated costs $\hat{c}$ as quadratic function.
For other convex versions of the approximated cost, a convex optimization with
the above linear constraints can be formulated.


Silvia Frassu$^{\sharp,*}$, Tongxing Li$^{\natural}$ Giuseppe Viglialoro$^{\sharp}$

$^\sharp$Dipartimento di Matematica e Informatica
Università di Cagliari
Via Ospedale 72, 09124. Cagliari (Italy)

$^{\natural}$School of Control Science and Engineering
Shandong University
Jinan, Shandong, 250061 (P. R. China)

We enter the details of two recent articles concerning as many chemotaxis models, one nonlinear and the other linear, and both with produced chemoattractant and saturated  chemorepellent.  More precisely, we are referring respectively to the papers “Boundedness in a nonlinear attraction-repulsion Keller–Segel system with production and consumption”, by S. Frassu, C. van der Mee and G. Viglialoro [J. Math. Anal. Appl. 504(2):125428, 2021] and “Boundedness in a chemotaxis system with consumed chemoattractant and produced chemorepellent”, by S. Frassu and G. Viglialoro [Nonlinear Anal. 213:112505, 2021]. These works, when properly analyzed, leave open room for some improvement of their results. We generalize the outcomes of the mentioned articles, establish other statements and put all the claims together; in particular, we select the sharpest ones and schematize them. Moreover, we complement our research also when logistic sources are considered in the overall study.

§ PREAMBLE

For details and discussions on the meaning of the forthcoming model, especially in the frame of chemotaxis phenomena and related variants, as well as for mathematical motivations and connected state of the art, we refer to [1, 2]. These articles will be often cited throughout this work.

§ PRESENTATION OF THE THEOREMS

Let $\Omega \subset \R^n$, $n \geq 2$, be a bounded and smooth domain, $\chi,\xi,\delta>0$, $m_1,m_2,m_3\in\R$, $f(u), g(u)$ and $h(u)$ be reasonably regular functions generalizing the prototypes $f(u)=K u^\alpha$, $g(u)=\gamma u^l$, and $h(u)=k u - \mu u^{\beta}$ with $K,\gamma, \mu>0$, $k \in \R$ and suitable $\alpha, l, \beta>0$. Once nonnegative initial configurations $u_0$ and $v_0$ are fixed, we aim at deriving sufficient conditions involving the above data so to ensure that the following attraction-repulsion chemotaxis model

\begin{equation}\label{problem}
\begin{cases}
u_t= \nabla \cdot ((u+1)^{m_1-1}\nabla u - \chi u(u+1)^{m_2-1}\nabla v+\xi u(u+1)^{m_3-1}\nabla w) + h(u) & \text{ in } \Omega \times (0,T_{\rm{max}}),\\
v_t=\Delta v-f(u)v  & \text{ in } \Omega \times (0,T_{\rm{max}}),\\
0= \Delta w - \delta w + g(u)& \text{ in } \Omega \times (0,T_{\rm{max}}),\\
u_{\nu}=v_{\nu}=w_{\nu}=0 & \text{ on } \partial \Omega \times (0,T_{\rm{max}}),\\
u(x,0)=u_0(x), \; v(x,0)=v_0(x) & x \in \bar\Omega,
\end{cases}
\end{equation}

admits classical solutions which are global and uniformly bounded in time. Specifically, we look for nonnegative functions $u=u(x,t), v=v(x,t), w=w(x,t)$ defined for $(x,t) \in \bar{\Omega}\times [0,T_{\rm{max}})$, and $T_{\rm{max}}=\infty$, with the properties that

\begin{equation}\label{ClassicalAndGlobability}
\begin{cases}
u,v\in C^0(\bar{\Omega}\times [0,\infty))\cap  C^{2,1}(\bar{\Omega}\times (0,\infty)), w\in C^0(\bar{\Omega}\times [0,\infty))\cap  C^{2,0}(\bar{\Omega}\times (0,\infty)),\\ (u,v,w) \in (L^\infty((0, \infty);L^{\infty}(\Omega)))^3,
\end{cases}
\end{equation}

and pointwisely satisfying all the relations in problem (<ref>).

To this scope, let  $f$, $g$ and $h$ be such that

\begin{equation}\label{f}
f,g \in C^1(\R) \quad \textrm{with} \quad   0\leq f(s)\leq Ks^{\alpha}  \textrm{ and } \gamma s^l\leq g(s)\leq \gamma s(s+1)^{l-1},\;  \textrm{for some}\; K,\,\gamma,\,\alpha>0,\, l\geq 1 \quad \textrm{and all } s \geq 0,
\end{equation}

\begin{equation}\label{h}
h \in C^1(\R) \quad \textrm{with} \quad h(0)\geq 0  \textrm{ and } h(s)\leq k s-\mu s^{\beta}, \quad  \textrm{for some}\quad k \in \R,\,\mu>0,\, \beta>1\, \quad \textrm{and all } s \geq 0.
\end{equation}

Then we prove these two theorems.

Let $\Omega$ be a smooth and bounded domain of $\mathbb{R}^n$, with $n\geq 2$, $\chi, \xi, \delta$ positive reals, $l \geq 1$ and $h \equiv 0$. Moreover, for $\alpha >0$ and $m_1, m_2, m_3 \in \R$, let $f$ and $g$ fulfill (<ref>) for each of the following cases:

*  $\alpha \in \left(0, \frac{1}{n}\right]$ and $m_1>\min\left\{2m_2-(m_3+l),\max\left\{2m_2-1,\frac{n-2}{n}\right\}, m_2 - \frac{1}{n}\right\}=:\mathcal{A}$,

*  $\alpha \in \left(\frac{1}{n},\frac{2}{n}\right)$ and $m_1>m_2 + \alpha - \frac{2}{n}=:\mathcal{B}$,

*  $\alpha \in \left[\frac{2}{n},1\right]$ and $m_1>m_2 + \frac{n\alpha-2}{n\alpha-1}=:\mathcal{C}$.

Then for any initial data $(u_0,v_0)\in (W^{1,\infty}(\Omega))^2$, with $u_0, v_0\geq 0$ on $\bar{\Omega}$, there exists a unique triplet $(u,v,w)$ of nonnegative functions, uniformly bounded in time and classically solving problem (<ref>).

Under the same hypotheses of Theorem <ref> and $\beta>1$, let $h$ comply with (<ref>).
Moreover, for $\alpha >0$ and $m_1, m_2, m_3 \in \R$, let $f$ and $g$ fulfill (<ref>) for each of the following cases:

*  $\alpha \in \left(0, \frac{1}{n}\right]$ and $m_1>\min\left\{2m_2-(m_3+l), 2m_2-\beta\right\}=:\mathcal{D}$,

*  $\alpha \in \left(\frac{1}{n},1\right)$ and $m_1>\min\left\{2m_2+1-(m_3+l), 2m_2+1-\beta\right\}=:\mathcal{E}$.

Then the same claim holds true.

When the logistic term $h$ does not take part in the model, problem (<ref>) has been already analyzed in [1] for the nonlinear diffusion and sensitivities case, and in [2] for the linear scenario; nevertheless, in these papers only small values of $\alpha$ are considered. Precisely, for $\alpha$ belonging to $(0,\frac{1}{2}+\frac{1}{n}),$  boundedness is ensured:

* in <cit.> for $m_1,m_2,m_3\in \R$ and $l\geq 1$, under the assumption

* in <cit.> for either $m_1=m_2=m_3=l=1$, under the assumption

\frac{2^\frac{2}{n}\frac{2}{n}^{n+1}(\frac{2}{n}-1)(n^2+n)}{(\frac{2}{n}+1)^{\frac{2}{n}+1}}\right)^\frac{2}{n}\lVert v_0\rVert_{L^\infty(\Omega)}^\frac{4}{n}=:\mathcal{G},$$

or in <cit.> for $m_1=m_2=m_3=1$ and any $l>1$.

In light of Theorems <ref> and <ref>, herein we develop an analysis dealing also with values of $\alpha$ larger than $\frac{1}{2}+\frac{1}{n}$. Additionally, for $\alpha$ belonging to some sub-intervals of $(0,\frac{1}{2}+\frac{1}{n})$ we improve <cit.> and  <cit.>. On the other hand, the introduction of $h$ allows us to obtain further generalizations and/or claims.

All this aspects are put together into Table <ref>. It, when possible, also indicates which of the assumptions taken from <cit.>, <cit.>, and Theorems <ref> and <ref> are the mildest leading to boundedness.

$m_2$    $m_3$   $l$   1c|$\alpha$   1c|$m_1$   $\chi$   $\xi$    Reference         Implication

$1$    $1$   $1$   1c|$[\frac{2}{n},1)$   1c|$1$   $\R^+$   1c|$ >\mathcal{G}$    Remark <ref>, generalizing <cit.>

$1$    $1$   $>1$   1c|$ [\frac{2}{n},1)$   1c|$1$   $ \R^+$   1c|$\R^+$    Remark <ref>, generalizing <cit.>

$ \R$    $\R$   $\geq 1$   1c|$ (\frac{1}{n},1)$   1c|$>\mathcal{F}$   $\R^+$   1c|$\R^+$    Remark <ref>, generalizing <cit.>

$ \R$    $\R$   $\geq 1$   1c|$ (0,\frac{1}{n}]$   1c|$>\mathcal{A}$   $\R^+$   1c|$\R^+$   Th. <ref>        $\ast$ and $\ast\ast$

$\R$     $\R$   $\geq 1$   1c|$(\frac{1}{n},\frac{2}{n})$   1c|$>\mathcal{B}$   $\R^+$   1c|$\R^+$   Th. <ref>        $\ast$

$\R$     $\R$   $\geq 1$   1c|$[\frac{2}{n},1]$   1c|$>\mathcal{C}$   $\R^+$   1c|$\R^+$   Th. <ref>

$m_2$    $m_3$   $l$   $\beta$   1c|$\alpha$   1c|$m_1$   $\chi$   $\xi$   $k$   $\mu$    Reference

$1$    $1$   $1$   $>2$   1c|$(\frac{1}{n},1)$   1c|$1$   $\R^+$   $\R^+$   $ \R$   1c|$\R^+$   Th. <ref>

$ \R$    $\R$   $\geq 1$   $>1$   1c|$ (0,\frac{1}{n}]$   1c|$>\mathcal{D}$   $\R^+$   $ \R^+$   $ \R$   1c|$\R^+$   Th. <ref>

$\R$     $\R$   $\geq 1$   $>1$   1c|$(\frac{1}{n},1)$   1c|$>\mathcal{E}$   $\R^+$   $ \R^+$   $ \R$   1c|$\R^+$   Th. <ref>

Schematization collecting the ranges of the parameters involved in model (<ref>) for which boundedness of its solutions is established for any fixed initial distribution $u_0$ and $v_0$. The symbol $\ast$ stands for “improves  <cit.> and recovers <cit.>”  and $\ast\ast$ for “improves <cit.>”. ($\mathcal{A}, \mathcal{B}, \mathcal{C}, \mathcal{D}, \mathcal{E}, \mathcal{F}$ are defined above.)

§ LOCAL WELL POSEDNESS, BOUNDEDNESS CRITERION, MAIN ESTIMATES AND ANALYSIS OF PARAMETERS

For $\Omega$, $\chi,\xi,\delta$, $m_1,m_2,m_3$ and $f, g, h$ as above, from now on with  $u, v, w \geq 0$ we refer to functions of $(x,t) \in \bar{\Omega}\times [0,T_{\rm{max}})$, for some finite $T_{\rm{max}}$, classically solving problem (<ref>) when nonnegative initial data $(u_0,v_0)\in (W^{1,\infty}(\Omega))^2$ are provided.  In particular, $u$ satisfies

\begin{equation}\label{massConservation}
\int_\Omega u(x, t)dx \leq m_0 \quad \textrm{for all }\, t \in (0,\TM),
\end{equation}

whilst $v$ is such that

\begin{equation*}\label{MaxPrincV}
0 \leq v\leq \lVert v_0\rVert_{L^\infty(\Omega)}\quad \textrm{in}\quad \Omega \times (0,T_{\rm{max}}).
\end{equation*}

Further, globality and boundedness of $(u,v,w)$ (in  the sense of (<ref>)) are ensured whenever (boundedness criterion) the $u$-component belongs to $L^\infty((0,T_{\rm{max}});L^p(\Omega))$, with $p>1$ arbitrarily large, and uniformly with respect $t\in (0,T_{\rm{max}})$.

These basic statements can be proved by standard reasoning; in particular, when $h\equiv 0$ they verbatim follow from <cit.> and relation (<ref>) is the well-known mass conservation property. Conversely, in the presence of the logistic terms $h$ as in (<ref>), some straightforward adjustments have to be considered and the $L^1$-bound of $u$ is consequence of an integration of the first equation in (<ref>) and an application of the Hölder inequality: precisely for $k_+=\max\{k,0\}$

\[
\frac{d}{dt} \int_\Omega u = \int_\Omega h(u) =k \int_\Omega u - \mu \int_\Omega u^{\beta} \leq k_+ \int_\Omega u - \frac{\mu}{|\Omega|^{\beta-1}} \left(\int_\Omega u\right)^{\beta}
\quad \textrm{for all }\, t \in (0,\TM),
\]

and we can conclude by invoking an ODI-comparison argument.

In our computations, beyond the above positions, some uniform bounds of $\|v(\cdot,t)\|_{W^{1,s}(\Omega)}$ are required. In this sense, the following lemma gets
the most out from $L^p$-$L^q$ (parabolic) maximal regularity; this is a cornerstone and for some small values of $\alpha$ the succeeding $W^{1,s}$-estimates are sharper than the $W^{1,2}$-estimates derived in [1, 2], and therein employed.

There exists $c_0>0$ such that $v$ fulfills

\begin{equation}\label{Cg}
\int_\Omega |\nabla v(\cdot, t)|^s\leq c_0 \quad \textrm{on } \,  (0,\TM)
\begin{cases}
\; \textrm{for all } s \in [1,\infty) & \textrm{if } \alpha \in \left(0, \frac{1}{n}\right],\\
\;  \textrm{for all } s \in \left[1, \frac{n}{(n\alpha-1)}\right) & \textrm{if } \alpha \in \left(\frac{1}{n},1\right].
\end{cases}
\end{equation}

For each $\alpha \in (0,1]$, there is $\rho >\frac{1}{2}$ such that for all $s \in \left[\frac{1}{\alpha},\frac{n}{(n\alpha-1)_+}\right)$ we have
$\frac{1}{2}<\rho <1-\frac{n}{2}\big(\alpha-\frac{1}{s}\big)$. From $1-\rho-\frac{n}{2}\big(\alpha-\frac{1}{s}\big)>0$, the claim follows invoking properties related to the Neumann heat semigroup, exactly as done in the second part of <cit.>.

We will make use of this technical result.

Let $n\in \N$, with $n\geq 2$, $m_1>\frac{n-2}{n}$, $m_2,m_3\in \R$  and $\alpha \in (0,1]$. Then there is $s \in [1,\infty)$, such that for proper
$p,q\in[1,\infty)$, $\theta$ and $\theta'$, $\mu$ and $\mu'$ conjugate exponents, we have that

\begin{align}
a_1&= \frac{\frac{m_1+p-1}{2}\left(1-\frac{1}{(p+2m_2-m_1-1)\theta}\right)}{\frac{m_1+p-1}{2}+\frac{1}{n}-\frac{1}{2}},
&  a_2&=\frac{q\left(\frac{1}{s}-\frac{1}{2\theta'}\right)}{\frac{q}{s}+\frac{1}{n}-\frac{1}{2}},  \nonumber \\
a_3 &= \frac{\frac{m_1+p-1}{2}\left(1-\frac{1}{2\alpha\mu}\right)}{\frac{m_1+p-1}{2}+\frac{1}{n}-\frac{1}{2}},
&  a_4  &= \frac{q\left(\frac{1}{s}-\frac{1}{2(q-1)\mu'}\right)}{\frac{q}{s}+\frac{1}{n}-\frac{1}{2}},\nonumber\\
\kappa_1 &  =\frac{\frac{p}{2}\left(1- \frac{1}{p}\right)}{\frac{m_1+p-1}{2}+\frac{1}{n}-\frac{1}{2}} \nonumber, & \kappa_2 & =  \frac{q - \frac{1}{2}}{q+\frac{1}{n}-\frac{1}{2}},
\end{align}

belong to the interval $(0,1)$. If, additionally,

\begin{equation}\label{Restrizionem1-m2-alphaPiccolo}
\alpha \in \left(0,\frac{1}{n}\right] \; \text{and}\; m_2-m_1<\frac{1}{n},
\end{equation}

\begin{equation}\label{Restrizionem1-m2-alphaGrande}
\alpha \in \left(\frac{1}{n},\frac{2}{n}\right) \; \text{and}\; m_2-m_1<\frac{2}{n}-\alpha,
\end{equation}

\begin{equation}\label{Restrizionem1-m2-alphaGrandeBis}
\alpha \in \left[\frac{2}{n},1\right] \; \text{and}\; m_2-m_1<\frac{2-n\alpha}{n\alpha-1},
\end{equation}

these futher relations hold true:

\begin{equation*}\label{MainInequalityExponents}
\beta_1 + \gamma_1 =\frac{p+2m_2-m_1-1}{m_1+p-1}a_1+\frac{1}{q}a_2\in (0,1) \;\textrm{ and }\;	\beta_2 + \gamma_2= \frac{2 \alpha }{m_1+p-1}a_3+\frac{q-1}{q}a_4
\in (0,1).
\end{equation*}

For any $s\geq 1$, let $\theta'>\max\left\{\frac{n}{2},\frac{s}{2}\right\}$  and $\mu>\max\left\{\frac{1}{2\alpha},\frac{n}{2}\right\}$. Thereafter, for

\begin{equation}\label{Prt_q}
\begin{cases}
q > \max \left\{\frac{(n-2)}{n}\theta', \frac{s}{2\mu'}+1\right\} \\
p>\max \left\{2-\frac{2}{n}-m_1,\frac{1}{\theta}-2m_2+m_1+1, \frac{(2m_2-m_1-1)(n-2)\theta-nm_1+n}{n-(n-2)\theta},\frac{2\alpha \mu(n-2)}{n} -m_1+1\right\},
\end{cases}
\end{equation}

it can be seen that $a_i,k_2\in (0,1)$, for any $i=1,2,3,4.$ On the other hand, $k_1\in (0,1)$ also thanks to the assumption $m_1>\frac{n-2}{n}.$

As to the second part, we distinguish three cases: $\alpha \in \left(0,\frac{1}{n}\right]$, $\alpha \in \left(\frac{1}{n},\frac{2}{n}\right)$ and
$\alpha \in \left[\frac{2}{n},1\right].$ (We insert Figure <ref> to clarify the proof, by focusing on the relation involving the values of $\alpha$, $s$ and
$\theta'$ in terms of assumptions (<ref>), (<ref>), (<ref>).)

$\circ$ $\alpha \in \left(0,\frac{1}{n}\right]$.
For $s>\frac{2\mu'}{2\mu'-1}$ arbitrarily large, consistently with (<ref>), we take $p=q=s$ and $\theta'=s\omega$, for some $\omega>\frac{1}{2}$. Computations provide

\begin{equation*}
0<	\beta_1+\gamma_1=\frac{s+2m_2-m_1-1-\frac{1}{\theta}}{m_1+s-2+\frac{2}{n}}+\frac{2-\frac{1}{\omega}}{s+\frac{2s}{n}},
\end{equation*}

\begin{equation*}
0<	\beta_2+\gamma_2=\frac{2\alpha-\frac{1}{\mu}}{m_1+s-2+\frac{2}{n}}+\frac{2s-2-\frac{s}{\mu'}}{s+\frac{2s}{n}}.
\end{equation*}

In light of the above positions, the largeness of $s$ infers $\theta$ arbitrarily close to $1$. Further, by choosing $\omega$ approaching $\frac{1}{2}$, continuity arguments imply that  $\beta_1+\gamma_1<1$ whenever restriction (<ref>) is satisfied, whereas $\beta_2+\gamma_2<1$ comes from $\mu>\frac{n}{2}.$

$\circ$ $\alpha \in \left(\frac{1}{n},\frac{2}{n}\right).$ First let $s$ be arbitrarily close to $\frac{n}{n\alpha-1}$ and let $q=\frac{p}{2}$ fulfill (<ref>). Then, in these circumstances it holds that $\max\left\{\frac{s}{2},\frac{n}{2}\right\}=\frac{s}{2}$, so that restriction on $\theta'$ reads $\theta'>\frac{s}{2}$.  Subsequently,

\begin{equation*}
0<	\beta_1+\gamma_1=\frac{p+2m_2-m_1-1-\frac{1}{\theta}}{m_1+p-2+\frac{2}{n}}+\frac{2-\frac{s}{\theta'}}{p+\frac{2s}{n}-s},
\end{equation*}

\begin{equation*}
0<	\beta_2+\gamma_2=\frac{2\alpha-\frac{1}{\mu}}{m_1+p-2+\frac{2}{n}}+\frac{p-2-\frac{s}{\mu'}}{p+\frac{2s}{n}-s}.
\end{equation*}

Since from $\theta'>\frac{s}{2}$ we have that $\theta'$ approaches $\frac{n}{2(n\alpha-1)}$, similar arguments used above imply that upon enlarging $p$ one can see that condition (<ref>) leads to  $\beta_1+\gamma_1<1$. On the other hand,  in order to have $\beta_2+\gamma_2<1$ we have to invoke the above constrain on $\mu$, i.e., $\mu>\frac{1}{2\alpha}.$

$\circ$ $\alpha \in \left[\frac{2}{n},1\right].$ We only have to consider in the previous case that now $\theta'>\frac{n}{2}$, so concluding thanks to

declare function=
func2()= (<= 1/2) * (1/2)   +
and(> 1/2, <= 2/2) * (-+2/2)     +
(>= 2/2) * (-(2*-2)/(2*-1))
func3()= (<= 1/3) * (1/3)   +
and(> 1/3, <= 2/3) * (-+2/3)     +
(>= 2/3) * (-(3*-2)/(3*-1))
func4()= (<= 1/4) * (1/4)   +
and(> 1/4, <= 2/4) * (-+2/4)     +
(>= 2/4) * (-(4*-2)/(4*-1))
func5()= (<= 1/5) * (1/5)   +
and(> 1/5, <= 2/5) * (-+2/5)     +
(>= 2/5) * (-(5*-2)/(5*-1))
func10()= (<= 1/10) * (1/10)   +
and(> 1/10, <= 2/10) * (-+2/10)     +
(>= 2/10) * (-(10*-2)/(10*-1))

axis x line=middle, axis y line=middle,
ymin=-1.5, ymax=0.6, ytick=-1,-0.8,...,0.6, ylabel=$m_2-m_1$,
xmin=0, xmax=1, xtick=0,0.2,...,1,
axis y line=center, axis x line=center,
width=12cm, height=12cm,axis on top
[thick, red] func2(x);
[thick, blue] func3(x);
[thick, cyan] func4(x);

[thick, green] func10(x);
(42,92)  node  $n=10$;
(83,100)  node  $n=4$;
(87,121)  node  $n=3$;
(83,178)  node  $n=2$;
(94,152)  node  $\alpha$;

[(-],  thick, black ] (0.1,8) – (20,8);
[(-),  thick, black] (20,8) – (40,8);
[[-],   thick, black ] (40,8) – (100,8);

(1.4,2.3) node   $\tiny{0}$;
(40,2.3) node   $\tiny{2/n}$;
(20,2.3) node   $\tiny{1/n}$;
(98,2.3) node   $\tiny{1}$;

(10,35) node$\tiny{s\in [1,\infty)}$;
(10,25) node$\tiny{\downarrow}$;
(10,17) node   $\tiny{\theta' \approx \infty}$;
(30,35) node  $\tiny{s\in [1,\frac{n}{n\alpha-1})}$;
(30,25) node$\tiny{\downarrow}$;
(30,17) node   $\tiny{\theta' \gtrapprox \frac{n}{2(n\alpha-1)}}$;
(70,35) node$\tiny{s\in [1,\frac{n}{n\alpha-1})}$;
(70,25) node$\tiny{\downarrow}$;
(70,17) node$\tiny{\theta' \gtrapprox n/2}$;

The colored lines, functions of $\alpha$, represent the supremum of the difference $m_2-m_1$ for some space dimension $n$. Moreover, for the sub-intervals $(0,1/n]$,  $(1/n,2/n)$ and $[2/n,1]$ of $\alpha$, the corresponding range of $s$ and choice of $\theta'$ are also indicated.

In view of its importance in the computations, we have to point out that from the above lemma $s$ can be chosen arbitrarily large only when $\alpha \in \left(0,\frac{1}{n}\right]$.
In particular, as we will see, in such an interval the terms  $\int_{\Omega} (u+1)^{p+2m_2-m_1-1} \vert \nabla v\lvert^2$ and $\int_{\Omega} (u+1)^{2\alpha} \vert \nabla v\lvert^{2(q-1)}$, appearing in our reasoning, can be treated in two alternative ways: either invoking the Young inequality or the Gagliardo–Nirenberg one.

§ A PRIORI ESTIMATES AND PROOF OF THE THEOREMS

§.§ The non-logistic case

Recalling the globality criterion mentioned in $\S$<ref>, let us define the functional $y(t):=\int_\Omega (u+1)^p + \int_\Omega |\nabla v|^{2q}$,
with $p,q>1$ properly large (and, when required, with $p=q$), and let us dedicate to derive the desired uniform bound of $\int_\Omega u^p$.

In the spirit of Remark <ref>, let us start by analyzing the evolution in time of the functional $y(t)$ by relying on the Young inequality.

Let $\alpha \in \left(0,\frac{1}{n}\right]$. If $m_1,m_2,m_3 \in \R$ comply with either $m_1>2m_2-(m_3+l)$ or $m_1>\max\left\{2m_2-1,\frac{n-2}{n}\right\}$,
then there exist $p, q>1$ such that $(u,v,w)$ satisfies for some $c_{16}, c_{17}, c_{18}>0$

\begin{equation}\label{MainInequality}
\frac{d}{dt} \left(\int_\Omega (u+1)^p + \int_\Omega  |\nabla v|^{2q}\right) + c_{16} \int_\Omega |\nabla |\nabla v|^q|^2 + c_{17} \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2 \leq c_{18} \quad \text{on } (0,\TM).
\end{equation}

Let $p=q>1$ sufficiently large; moreover, in view of Remark <ref>, from now on, when necessary we will tacitly enlarge these parameters.

In the first part of the proof we focus on the estimate of the term $\frac{d}{dt} \int_\Omega (u+1)^p$.
Standard testing procedures provide

\begin{equation*}
\begin{split}
\frac{d}{dt} \int_\Omega (u+1)^p =\int_\Omega p(u+1)^{p-1}u_t &= -p(p-1) \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2
+p(p-1)\chi \int_\Omega u(u+1)^{m_2+p-3} \nabla u \cdot \nabla v \\
&-p(p-1)\xi \int_\Omega u(u+1)^{m_3+p-3} \nabla u \cdot \nabla w \quad \text{on } (0,\TM).
\end{split}
\end{equation*}

By reasoning as in <cit.>, we obtain for $\epsilon_1, \epsilon_2, \tilde{\sigma}$ positive, and for all $t \in (0,\TM)$, some $c_1>0$ such that

\begin{equation}\label{Estim_1_For_u^p}
\begin{split}
\frac{d}{dt} \int_\Omega (u+1)^p & \leq -p(p-1) \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 +p(p-1)\chi \int_\Omega u(u+1)^{m_2+p-3} \nabla u \cdot \nabla v \\
&+ \left(\epsilon_1 + \tilde{\sigma} -  \frac{p(p-1)\xi \gamma}{2^{p-1}(m_3+p-1)} \right) \int_\Omega (u+1)^{m_3+p+l-1}
+ \epsilon_2 \int_\Omega |\nabla (u+1)^\frac{m_1+p-1}{2}|^2 + c_1.
\end{split}
\end{equation}

Let us now discuss the cases $m_1>2m_2-(m_3+l)$ and $m_1>\max\left\{2m_2-1,\frac{n-2}{n}\right\}$, respectively.
A double application of the Young inequality in (<ref>) and bound (<ref>) give

\begin{equation}\label{Young}
\begin{split}
&p(p-1)\chi \int_\Omega u (u+1)^{m_2+p-3} \nabla u \cdot \nabla v \leq \epsilon_3 \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2
+ c_2 \int_\Omega (u+1)^{p+2m_2-m_1-1} |\nabla v|^2\\
& \leq \epsilon_3 \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 + \epsilon_4 \int_\Omega |\nabla v|^{s} + c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}}\\
& \leq \epsilon_3 \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 + c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}} + c_4 \quad \text{on }\, (0,\TM),
\end{split}
\end{equation}

with $\epsilon_3, \epsilon_4 >0$ and some positive $c_2, c_3, c_4$.
From $m_1> 2m_2-(m_3+l)$, we have $\frac{(p+2m_2-m_1-1)s}{s-2} < (m_3+p+l-1)$, and for every $\epsilon_5>0$, Young's inequality yields some $c_5>0$ entailing

\begin{equation}\label{Young1}
c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}} \leq \epsilon_5 \int_\Omega (u+1)^{m_3+p+l-1} + c_5 \quad \text{for all }\, (0,\TM).
\end{equation}

Now, we note that $m_1>2 m_2-1$ implies $\frac{(p+2m_2-m_1-1)s}{s-2} < p$, and the Young inequality allows us to rephrase (<ref>) in an alternative way:

\begin{equation}\label{Young4}
c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}} \leq \epsilon_6 \int_\Omega (u+1)^p + c_6 \quad \text{on }\, (0,\TM),
\end{equation}

with $\epsilon_6>0$ and positive $c_6$.
Further, an application of the Gagliardo–Nirenberg inequality and property (<ref>) yield

\begin{equation*}\label{Theta_2}
\theta=\frac{\frac{n(m_1+p-1)}{2}\left(1-\frac{1}{p}\right)}{1-\frac{n}{2}+\frac{n(m_1+p-1)}{2}}\in (0,1),
\end{equation*}

so giving for $c_7, c_8>0$

\begin{equation*}
\begin{split}
\int_{\Omega} (u+1)^p&= \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2p}{m_1+p-1}}(\Omega)}^{\frac{2p}{m_1+p-1}}\\
&\leq c_7 \|\nabla (u+1)^{\frac{m_1+p-1}{2}}\|_{L^2(\Omega)}^{\frac{2p}{m_1+p-1}\theta} \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{2p}{m_1+p-1}(1-\theta)} + c_7 \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{2p}{m_1+p-1}}\\
& \leq c_8 \Big(\int_\Omega |\nabla (u+1)^\frac{m_1+p-1}{2}|^2\Big)^{\kappa_1}+ c_8 \quad \text{ for all } t \in(0,\TM).
\end{split}
\end{equation*}

Since $\kappa_1 \in (0,1)$ (see Lemma <ref>), for any positive $\epsilon_7$  thanks to the Young inequality we arrive for some positive $c_9>0$ at

\begin{equation}\label{GN2}
\begin{split}
\epsilon_6 \int_\Omega (u+1)^p  &\leq \epsilon_7 \int_\Omega |\nabla (u+1)^\frac{m_1+p-1}{2}|^2  + c_9 \quad \text{on } (0,\TM).
\end{split}
\end{equation}

By plugging estimate (<ref>) into relation (<ref>), and by relying on bound (<ref>) (or, alternatively to (<ref>), relations (<ref>) and (<ref>)), infer for appropriate $\tilde{\epsilon}_1, \tilde{\epsilon}_2>0$ and proper $c_{10}>0$

\begin{equation}\label{ClaimU}
\begin{split}
\frac{d}{dt} \int_\Omega (u+1)^p &\leq \left(-\frac{4p(p-1)}{(m_1+p-1)^2} + \tilde{\epsilon}_1 \right) \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\\
&+ \left(\tilde{\epsilon}_2 - \frac{p(p-1)\xi \gamma}{2^{p-1}(m_3+p-1)}\right) \int_\Omega (u+1)^{m_3+p+l-1} + c_{10} \quad \text{for all } t \in (0,\TM),
\end{split}
\end{equation}

where we also exploited that

\begin{equation}\label{GradU}
\int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 = \frac{4}{(m_1+p-1)^2} \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2 \quad  \text{on } (0,\TM).
\end{equation}

Now, as to the term $\frac{d}{dt} \int_\Omega  |\nabla v|^{2q}$ of the functional $y(t)$, reasoning similarly as in <cit.>, we obtain for some
$c_{11}, c_{12}>0$

\begin{equation}\label{Estim_gradV}
\frac{d}{dt}\int_\Omega  |\nabla v|^{2q}+ q \int_\Omega |\nabla v|^{2q-2} |D^2v|^2 \leq c_{11} \int_\Omega u^{2\alpha} |\nabla v|^{2q-2} +c_{12} \quad \textrm{on } \;  (0,\TM).
\end{equation}

Moreover, Young's inequalities and bound (<ref>) give for every arbitrary $\epsilon_8, \epsilon_9>0$ and some $c_{13}, c_{14}, c_{15}>0$

\begin{equation}\label{Estimat_nablav^2q+2}
\begin{split}
&c_{11} \int_\Omega u^{2\alpha} |\nabla v|^{2q-2} \leq \epsilon_8 \int_\Omega (u+1)^{m_3+p+l-1} + c_{13} \int_\Omega |\nabla v|^{\frac{2(q-1)(m_3+p+l-1)}{m_3+p+l-1-2\alpha}}\\
& \leq \epsilon_8 \int_\Omega (u+1)^{m_3+p+l-1} + \epsilon_9 \int_\Omega |\nabla v|^s + c_{14}
\leq \epsilon_8 \int_\Omega (u+1)^{m_3+p+l-1} + c_{15} \quad \textrm{for all} \quad t \in (0,\TM).
\end{split}
\end{equation}

Therefore, by inserting relation (<ref>) into (<ref>) and adding (<ref>), we have the claim for a proper choice of
$\tilde{\epsilon}_1, \tilde{\epsilon}_2, \epsilon_8$ and some positive $c_{16}, c_{17}, c_{18}$, also by taking into account the relation (see <cit.>)

\begin{equation}\label{GradV}
\vert \nabla \lvert \nabla v\rvert^q\rvert^2=\frac{q^2}{4}\lvert \nabla v \rvert^{2q-4}\vert \nabla \lvert \nabla v\rvert^2\rvert^2=q^2\lvert \nabla v \rvert^{2q-4}\lvert D^2v \nabla v \rvert^2\leq q^2|\nabla v|^{2q-2} |D^2v|^2.
\end{equation}

Let us now turn our attention when, as mentioned before, the Gagliardo–Nirenberg inequality is employed. In this case, we can derive information not only for
$\alpha \in \left(0,\frac{1}{n}\right]$ but also for $\alpha \in \left(\frac{1}{n},1\right]$.

If $m_1,m_2\in \R$ and $\alpha>0$ are taken accordingly to (<ref>), (<ref>),
(<ref>), then there exist $p, q>1$ such that $(u,v,w)$ satisfies a similar inequality as in (<ref>).

For $s$, $p$ and $q$ taken accordingly to Lemma <ref> (in particular, $p=q$ for $\alpha \in \left(0,\frac{1}{n}\right]$, and $q=\frac{p}{2}$
for $\alpha \in \left(\frac{1}{n},1\right]$), let $\theta, \theta', \mu, \mu'$, $a_1,a_2, a_3, a_4$ and $\beta_1, \beta_2, \gamma_1, \gamma_2$ be therein defined.

With a view to Lemma <ref>, by manipulating relation (<ref>) and focusing on the first inequality in (<ref>) and on (<ref>), proper $\epsilon_1, \tilde{\sigma}$ lead to

\begin{equation}\label{Somma}
\begin{split}
&\frac{d}{dt} \left(\int_\Omega (u+1)^p + \int_\Omega  |\nabla v|^{2q}\right) + q \int_\Omega |\nabla v|^{2q-2} |D^2v|^2
\leq \left(-\frac{4p(p-1)}{(m_1+p-1)^2} + \tilde{\epsilon}_1 \right) \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\\
&+  c_2 \int_\Omega (u+1)^{p+2m_2-m_1-1} |\nabla v|^2 + c_{11} \int_\Omega u^{2\alpha} |\nabla v|^{2q-2} +c_{19} \quad \textrm{on } \;  (0,\TM),
\end{split}
\end{equation}

for some $c_{19}>0$ (we also used relation (<ref>)).
In this way, we can estimate the second and third integral on the right-hand side of (<ref>) by applying the Hölder inequality so to have

\begin{equation} \label{H1}
\int_{\Omega} (u+1)^{p+2m_2-m_1-1} |\nabla v|^2 \leq  \left(\int_{\Omega} (u+1)^{(p+2m_2-m_1-1)\theta}\right)^{\frac{1}{\theta}}
\left(\int_{\Omega} |\nabla v|^{2 \theta'}\right)^{\frac{1}{\theta'}} \quad \textrm{ on } (0, \TM),
\end{equation}

\begin{equation} \label{H2}
\int_{\Omega} (u+1)^{2\alpha} |\nabla v|^{2q-2} \leq
\left(\int_{\Omega} (u+1)^{2\alpha\mu}\right)^{\frac{1}{\mu}} \left(\int_{\Omega} |\nabla v|^{2(q-1)\mu'}\right)^{\frac{1}{\mu'}}\quad  \textrm{ on }  (0,\TM).
\end{equation}

By invoking the Gagliardo–Nirenberg inequality and bound (<ref>), we obtain for some $c_{20}, c_{21}>0$

\begin{align}\label{a1}
&\left(\int_{\Omega} (u+1)^{(p+2m_2-m_1-1)\theta}\right)^{\frac{1}{\theta}}= \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2(p+2m_2-m_1-1)}{m_1+p-1}\theta}(\Omega)}^{\frac{2(p+2m_2-m_1-1)}{m_1+p-1}}\\ \nonumber
& \leq c_{20} \|\nabla(u+1)^{\frac{m_1+p-1}{2}}\|_{L^2(\Omega)}^{\frac{2(p+2m_2-m_1-1)}{m_1+p-1} a_1} \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{2(p+2m_2-m_1-1)}{m_1+p-1} (1-a_1)} + c_{20} \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{2(p+2m_2-m_1-1)}{m_1+p-1}} \\ \nonumber
&\leq c_{21} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_1}+ c_{21} \quad \textrm{ for all } \,t\in(0,\TM),
\end{align}

and for some $c_{22}, c_{23}>0$

\begin{align}\label{a3}
&\left(\int_{\Omega} (u+1)^{2\alpha\mu}\right)^{\frac{1}{\mu}}= \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{4\alpha\mu}{m_1+p-1}}(\Omega)}^{\frac{4\alpha}{m_1+p-1}}\\ \nonumber
& \leq c_{22} \|\nabla(u+1)^{\frac{m_1+p-1}{2}}\|_{L^2(\Omega)}^{\frac{4\alpha}{m_1+p-1}  a_3} \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{4\alpha}{m_1+p-1} (1-a_3)} + c_{22} \|(u+1)^{\frac{m_1+p-1}{2}}\|_{L^{\frac{2}{m_1+p-1}}(\Omega)}^{\frac{4\alpha}{m_1+p-1}}\\ \nonumber
&\leq c_{23} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_2}+ c_{23}
\quad \textrm{for all }\, t \in (0,\TM).
\end{align}

In a similar way, we can again apply the Gagliardo–Nirenberg inequality and bound (<ref>) and get for some $c_{24}, c_{25}>0$

\begin{align}\label{a2}
&\left(\int_{\Omega} |\nabla v|^{2 \theta'}\right)^{\frac{1}{\theta'}} =\| |\nabla v|^q \|_{L^{\frac{2 \theta'}{q}}(\Omega)}^{\frac{2}{q}}\leq c_{24} \|\nabla |\nabla v|^q\|_{L^2(\Omega)}^{\frac{2}{q}a_2} \| |\nabla v|^q\|_{L^{\frac{s}{q}}(\Omega)}^{\frac{2}{q}(1-a_2)} + c_{24} \| |\nabla v|^q\|_{L^{\frac{s}{q}}(\Omega)}^{\frac{2}{q}} \\ \nonumber
& \leq c_{25} \left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_1}+ c_{25}
\quad \textrm{ for all } t\in (0, \TM),
\end{align}

and for some $c_{26}, c_{27}>0$

\begin{align}\label{a4}
&\left(\int_{\Omega} |\nabla v|^{2(q-1) \mu'}\right)^{\frac{1}{\mu'}} =\| |\nabla v|^q \|_{L^{\frac{2(q-1)}{q}\mu'}(\Omega)}^{\frac{2(q-1)}{q}}
\leq c_{26} \|\nabla |\nabla v|^q\|_{L^2(\Omega)}^{\frac{2(q-1)}{q}  a_4} \| |\nabla v|^q\|_{L^{\frac{s}{q}}(\Omega)}^{\frac{2(q-1)}{q}  (1-a_4)}
+ c_{26} \| |\nabla v|^q\|_{L^{\frac{s}{q}}(\Omega)}^{\frac{2(q-1)}{q}} \\ \nonumber
&\leq c_{27} \left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_2}+ c_{27} \quad \textrm{ for every } t\in (0,\TM).
\end{align}

By plugging (<ref>) and (<ref>) into (<ref>) and taking into account (<ref>), (<ref>), (<ref>), (<ref>), we deduce for a proper
$\tilde{\epsilon}_1$, once inequality (<ref>) is considered, the following estimate for some $c_{28}, c_{29}, c_{30}, c_{31}, c_{32}>0$:

\begin{equation}\label{Somma1}
\begin{split}
&\frac{d}{dt} \left(\int_\Omega (u+1)^p + \int_\Omega  |\nabla v|^{2q}\right) + c_{28} \int_\Omega |\nabla |\nabla v|^q|^2
+ c_{29} \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2 - c_{30}\\
&\leq c_{31} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_1}
\left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_1} + c_{31} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_1}\\
& + c_{31} \left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_1}
+ c_{32} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_2} \left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_2}\\
&+ c_{32} \left(\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2\right)^{\beta_2} + c_{32} \left(\int_{\Omega} |\nabla |\nabla v|^q|^2 \right)^{\gamma_2} \quad \text{on }
\end{split}
\end{equation}

Since by Lemma <ref> we have that $\beta_1 + \gamma_1 <1$ and $\beta_2 + \gamma_2 <1$, and in particular $\beta_1, \gamma_1, \beta_2, \gamma_2 \in (0,1)$, we can treat the two integral products and the remaining four addenda of the right-hand side in a such way that eventually they are absorbed by the two integral terms involving the gradients in the left one. More exactly, to the products we apply

\[
a^{d_1}b^{d_2} \leq \epsilon(a+b)+c \quad \textrm{with } a,b\geq0, d_1,d_2 >0 \; \textrm{such that } d_1+d_2<1, \; \textrm{for all } \epsilon>0 \; \textrm{and some } c>0
\]

(achievable by means of applications of Young's inequality), and to the other terms the Young inequality. In this way, the resulting linear combination of $\int_{\Omega} |\nabla |\nabla v|^q|^2$ and $\int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2$ can be turned into $\frac{c_{28}}{2} \int_{\Omega} |\nabla |\nabla v|^q|^2 + \frac{c_{29}}{2} \int_{\Omega} |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2$, which coming back to (<ref>) infers what claimed.

§.§ The logistic case

For the logistic case we retrace part of the computations above connected to the usage of the Young inequality only.

If $m_1,m_2,m_3 \in \R$ comply with $m_1>2m_2-(m_3+l)$ or $m_1>2m_2-\beta$ whenever $\alpha \in \left(0, \frac{1}{n}\right]$, or
$m_1>2m_2 +1-(m_3+l)$ or $m_1>2m_2+1-\beta$ whenever $\alpha \in \left(\frac{1}{n},1\right)$, then there exist $p,q>1$ such that $(u,v,w)$ satisfies
a similar inequality as in (<ref>).

As in Lemma <ref>, in view of inequality (<ref>) and the properties of the logistic $h$ in (<ref>), relation (<ref>) now becomes for some positive $c_{33}$

\begin{equation}\label{Estim_1_For_u^p1}
\begin{split}
\frac{d}{dt} \int_\Omega (u+1)^p & \leq (-p(p-1)+\epsilon_3) \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 +c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}}\\
&+ \left(\epsilon_1 + \tilde{\sigma} -  \frac{p(p-1)\xi \gamma}{2^{p-1}(m_3+p-1)} \right) \int_\Omega (u+1)^{m_3+p+l-1}
+ \epsilon_2 \int_\Omega |\nabla (u+1)^\frac{m_1+p-1}{2}|^2\\
&+ pk_+ \int_\Omega (u+1)^p - p \mu  \int_\Omega (u+1)^{p-1}u^{\beta} + c_{33} \quad \text{ for all } t \in (0,\TM).
\end{split}
\end{equation}

Applying the inequality $(A+B)^p \leq 2^{p-1} (A^p+B^p)$ with $A,B \geq 0$ and $p>1$ to the last integral in (<ref>), implies that
$-u^{\beta} \leq -\frac{1}{2^{\beta-1}} (u+1)^{\beta}+1$; therefore

\begin{equation}\label{beta}
- p \mu  \int_\Omega (u+1)^{p-1} u^{\beta} \leq -\frac{p \mu}{2^{\beta-1}} \int_\Omega (u+1)^{p-1+\beta} + p \mu \int_\Omega (u+1)^{p-1} \quad \text{on } (0,\TM).
\end{equation}

Henceforth, by taking into account the Young inequality, we have that for $t \in (0,\TM)$

\begin{equation} \label{k}
pk_+ \int_\Omega (u+1)^p \leq \delta_1 \int_\Omega (u+1)^{p-1+\beta} + c_{34} \quad \text{and} \quad
p \mu \int_\Omega (u+1)^{p-1} \leq \delta_2 \int_\Omega (u+1)^{p-1+\beta} + c_{35},
\end{equation}

with $\delta_1, \delta_2 >0$ and some $c_{34}, c_{35} >0$.

Case $1$: $\alpha \in \left(0, \frac{1}{n}\right]$ and $m_1>2m_2-(m_3+l)$ or $m_1>2m_2-\beta$.
For $m_1> 2m_2-(m_3+l)$ we refer to Lemma <ref> and we take in mind inequality (<ref>).
Conversely, when $m_1>2 m_2-\beta$, we have that (recall $s$ may be arbitrary large) $\frac{(p+2m_2-m_1-1)s}{s-2} < p-1+\beta$, and by means of the Young inequality estimate (<ref>) can alternatively read

\begin{equation}\label{young4}
c_3 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)s}{s-2}} \leq \delta_3 \int_\Omega (u+1)^{p-1+\beta} + c_{36} \quad \text{on }\, (0,\TM),
\end{equation}

with $\delta_3>0$ and positive $c_{36}$.
By inserting estimates (<ref>) and (<ref>) into relation (<ref>), as well as taking into account (<ref>)
(or, alternatively to (<ref>), bound (<ref>)), for suitable $\hat{\epsilon}, \tilde{\epsilon}_2, \tilde{\delta} >0$ and some $c_{37}>0$ we arrive at

\begin{equation*}\label{ClaimUlog}
\begin{split}
\frac{d}{dt} \int_\Omega (u+1)^p &\leq \left(-\frac{4p(p-1)}{(m_1+p-1)^2} + \hat{\epsilon}\right) \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2
+ \left(\tilde{\epsilon_2} - \frac{p(p-1)\xi \gamma}{2^{p-1}(m_3+p-1)} \right) \int_\Omega (u+1)^{m_3+p+l-1}\\
& + \left(\tilde{\delta} - \frac{p \mu}{2^{\beta-1}} \right) \int_\Omega (u+1)^{p-1+\beta} + c_{37} \quad \text{ for all } t \in (0,\TM),
\end{split}
\end{equation*}

where we used again relation (<ref>). We conclude reasoning exactly as in the second part of the proof of Lemma <ref> and by choosing suitable
$\hat{\epsilon}, \tilde{\epsilon}_2, \tilde{\delta}, \epsilon_8$.

Case $2$: $\alpha \in \left(\frac{1}{n},1\right)$ and $m_1>2m_2 +1-(m_3+l)$ or $m_1>2m_2+1-\beta$.
Accordingly to Remark <ref>, since now $s$ cannot increase arbitrarily, relations (<ref>) and (<ref>) have to be differently discussed.
In particular, for some $\bar{c}_1>0$ we can estimate relation (<ref>) as follows:

\begin{equation*}
\begin{split}
&p(p-1)\chi \int_\Omega u (u+1)^{m_2+p-3} \nabla u \cdot \nabla v \leq \epsilon_3 \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2
+ c_2 \int_\Omega (u+1)^{p+2m_2-m_1-1} |\nabla v|^2\\
& \leq \epsilon_3 \int_\Omega (u+1)^{p+m_1-3} |\nabla u|^2 + \epsilon_4 \int_\Omega |\nabla v|^{2(p+1)} + \bar{c}_1 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)(p+1)}{p}}
\quad \text{on }\, (0,\TM).
\end{split}
\end{equation*}

Now, if $m_1>2m_2+1-(m_3+l)$, then some $p$ sufficiently large infers to $\frac{(p+2m_2-m_1-1)(p+1)}{p}< p+m_3+l-1$, so that for any positive $\bar{\epsilon}_1$
and some $\bar{c}_2>0$ we have

\[
\bar{c}_1 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)(p+1)}{p}} \leq \bar{\epsilon}_1 \int_\Omega (u+1)^{p+m_3+l-1} + \bar{c}_2 \quad \text{on } (0,\TM).
\]

Conversely, and in a similar way, for $m_1>2m_2+1-\beta$ we have for any positive $\bar{\epsilon}_2$ and some $\bar{c}_3>0$

\[
\bar{c}_1 \int_\Omega (u+1)^{\frac{(p+2m_2-m_1-1)(p+1)}{p}} \leq \bar{\epsilon}_2 \int_\Omega (u+1)^{p-1+\beta} + \bar{c}_3 \quad \text{on } (0,\TM).
\]

The remaining part of the proof follows as the previous case, by taking into account <cit.> for the term dealing with
$ \int_\Omega |\nabla v|^{2(p+1)}$.

As a by-product of what now obtained we are in a position to conclude.

§.§ Proof of Theorems <ref> and <ref>

Let $(u_0,v_0) \in (W^{1,\infty}(\Omega))^2$ with $u_0, v_0 \geq 0$ on $\bar{\Omega}$. For $f$ and $g$ as in (<ref>) and, respectively, for $f$, $g$ as in (<ref>) and
$h$ as in (<ref>), let $\alpha >0$ and let $m_1,m_2,m_3 \in \R$ comply with <ref>, <ref> and <ref>, respectively, <ref> and <ref>. Then, we refer to Lemmas <ref> and <ref>, respectively, Lemma <ref> and obtain for some

\begin{equation}\label{Estim_general_For_y_2}
y'(t) + C_1 \int_\Omega |\nabla (u+1)^{\frac{m_1+p-1}{2}}|^2 + C_2 \int_\Omega |\nabla |\nabla v|^q|^2 \leq C_3
\quad \text{ on } (0, \TM).
\end{equation}

Successively, the Gagliardo–Nirenberg inequality again makes that some positive constants $c_{38}, c_{39}$ imply from the one hand

\begin{equation*}\label{G_N2}
\int_\Omega (u+1)^p \leq c_{38} \Big(\int_\Omega |\nabla (u+1)^\frac{m_1+p-1}{2}|^2\Big)^{\kappa_1}+ c_{38}  \quad \text{for all } t \in (0,\TM),
\end{equation*}

(as already done in (<ref>)), and from the other

\begin{equation}\label{Estim_Nabla nabla v^q}
\int_\Omega \lvert \nabla v\rvert^{2q}=\lvert \lvert \lvert \nabla v\rvert^q\lvert \lvert_{L^2(\Omega)}^2
\leq c_{39} \lvert \lvert\nabla  \lvert \nabla v \rvert^q\rvert \lvert_{L^2(\Omega)}^{2\kappa_2} \lvert \lvert\lvert \nabla v \rvert^q\lvert \lvert_{L^\frac{1}{q}(\Omega)}^{2(1-\kappa_2)} +c_{39} \lvert \lvert \lvert \nabla v \rvert^q\lvert \lvert^2_{L^\frac{1}{q}(\Omega)}\quad \textrm{on } (0,\TM),
\end{equation}

with $\kappa_2$ already defined in Lemma <ref>.
Subsequently, the $L^s$-bound of $\nabla v$ in (<ref>) infers some $c_{40}>0$ such that

\begin{equation*}\label{Estim_Nabla nabla v^p^2}
\int_\Omega \lvert \nabla v\rvert^{2q}\leq c_{40} \Big(\int_\Omega \lvert \nabla \lvert \nabla v \rvert^q\rvert^2\Big)^{\kappa_2}+c_{40} \quad \text{for all } t \in (0,\TM).
\end{equation*}

In the same flavour of <cit.>, by using the estimates involving $\int_\Omega (u+1)^p$ and $\int_\Omega \rvert\nabla v\lvert^{2q}$,
relation (<ref>) provides positive constants  $c_{41}$ and $c_{42}$, and $\kappa=\min\{\frac{1}{\kappa_1},\frac{1}{\kappa_2}\}$ such that

\begin{equation*}\label{MainInitialProblemWithM}
\begin{cases}
y'(t)\leq c_{41}-c_{42} y^{\kappa}(t)\quad \textrm{for all } t \in (0,T_{\rm{max}}),\\
y(0)=\int_\Omega (u_0+1)^p+ \int_\Omega |\nabla v_0|^{2q}.
\end{cases}
\end{equation*}

Finally, ODE comparison principles imply $u \in L^\infty((0,\TM);L^p(\Omega))$, and the conclusion is a consequence of the boundedness criterion in

In the proofs of <cit.> and <cit.>, it is seen that the $L^2$ uniform estimate of $\nabla v$ is used to control some integral on $\partial \Omega$ (and this allows us to avoid to restrict our study to convex domains), as well as to deal with the term
$\int_\Omega \rvert \nabla v \lvert^{2p}$ with the Gagliardo–Nirenberg inequality; for instance we are referring to <cit.>, respectively.
Such finiteness of $\int_\Omega \rvert \nabla v \lvert^2$ is related to the values of $\alpha$ in these articles:
$\alpha \in \left(0, \frac{1}{2}+\frac{1}{n}\right)$ (see <cit.>). Apparently only $\nabla v \in L^\infty((0,\TM);L^1(\Omega))$ suffices to address these issues. Indeed, as far as the topological property of $\Omega$ is concerned, we can invoke <cit.> with $s=1$; on the other hand, for the question tied to the employment of the Gagliardo–Nirenberg inequality, we may operate as done in (<ref>). As a consequence, in view of Lemma <ref>, we have that  $\nabla v \in L^\infty((0,\TM);L^1(\Omega))$, so that <cit.> and <cit.> hold true for any $\alpha \in (0,1)$.

§.§.§ Acknowledgments

SF and GV are members of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM) and are partially supported by the research project Evolutive and Stationary Partial Differential Equations with a Focus on Biomathematics, funded by Fondazione di Sardegna (2019). GV is partially supported by MIUR (Italian Ministry of Education, University and Research) Prin 2017 Nonlinear Differential Problems via Variational, Topological and Set-valued Methods (Grant Number: 2017AYM8XW). TL is partially supported by NNSF of P. R. China (Grant No. 61503171), CPSF (Grant No. 2015M582091), and NSF of Shandong Province (Grant No. ZR2016JL021).

[1]
S. Frassu, C. van der Mee, and G. Viglialoro.
Boundedness in a nonlinear attraction-repulsion Keller–Segel
system with production and consumption.
J. Math. Anal. Appl., 504(2):125428, 2021.

[2]
S. Frassu and G. Viglialoro.
Boundedness in a chemotaxis system with consumed chemoattractant and
produced chemorepellent.
Nonlinear Anal., 213:112505, 2021.

[3]
S. Ishida, K. Seki, and T. Yokota.
Boundedness in quasilinear Keller–Segel systems of
parabolic-parabolic type on non-convex bounded domains.
J. Differential Equations, 256(8):2993–3010, 2014.