Bifurcation Analysis and Multiobjective Nonlinear Model Predictive Control of Drug Addiction Models

Lakshmi N Sridhar; Lakshmi N Sridhar

ISSN: 2455-2976

Journal of Cardiovascular Medicine and Cardiology

Review Article Open Access Peer-Reviewed

Bifurcation Analysis and Multiobjective Nonlinear Model Predictive Control of Drug Addiction Models

Lakshmi N Sridhar*

Chemical Engineering Department, University of Puerto Rico, Mayaguez, PR 00681, USA

Author and article information

*Corresponding author: Lakshmi N Sridhar, Chemical Engineering Department, University of Puerto Rico, Mayaguez, PR 00681, USA, E-mail: [email protected]

doi : 10.17352/2455-2976.000215

Received: 16 December, 2024 | Accepted: 23 December, 2024 | Published: 24 December, 2024

Keywords: Drug; Addiction; Bifurcation; Optimal control

Cite this as

Sridhar LN. Bifurcation Analysis and Multiobjective Nonlinear Model Predictive Control of Drug Addiction Models. J Cardiovasc Med Cardiol. 2024;11(4):096-102. Available from: 10.17352/2455-2976.000215

Copyright License

© 2024 Sridhar LN. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Bifurcation analysis and nonlinear model predictive control were performed on drug addiction models. Rigorous proof showing the existence of bifurcation (branch) points is presented along with computational validation. It is also demonstrated (both numerically and analytically) that the presence of the branch points was instrumental in obtaining the Utopia solution when the multiobjective nonlinear model prediction calculations were performed. Bifurcation analysis was performed using the MATLAB software MATCONT while the multi-objective nonlinear model predictive control was performed by using the optimization language PYOMO.

Main article text

Introduction

Mental health has become a significant focus for researchers and medical doctors in the last decade. Ironically, drug addiction is both cause and effect for the existence of mental health problems. People with mental health issues resort to drugs and drugs in turn lead to mental health problems. Additionally, drug addiction has led to a considerable amount of poverty and crime. It is therefore important to develop strategies to curb drug addiction. The problem of drug addiction has led to computational research to develop reliable techniques to be able to control drug addiction. This work aims to perform bifurcation analysis in conjunction with multiobjective nonlinear model predictive control (MNLMPC) calculations on models involving drug addiction. This paper is organized as follows. First, the background section with the literature review is presented. The bifurcation analysis techniques and the multiobjective nonlinear model predictive control strategies are presented followed by a description of how the presence of singular points affects the MNLMPC calculations. Two drug addiction example problems where MNLMPC calculations are performed in conjunction with bifurcation analysis are presented. It is numerically demonstrated that the presence of bifurcation points in the drug addiction models enables the MNLMPC calculations to converge to the Utopia solution. It is essential to develop rigorous strategies to be able to effectively combat the drug addiction problem minimizing the number of drug addicts and avoiding the wastage of resources. This paper is a step in that direction.

Background

Bae [1] studied the dynamics of tobacco addiction models. Mushayabasa, and co-workers [2-4] performed dynamic and optimal control studies of drug addiction models. Hasan, et al. [5] investigated the effect of having drug rehabilitation centers to combat drug addiction. Islam, et al. [6,7] developed a mathematical analysis of some dynamic Models of drug addiction, while Lavi, et al. [8] studied the dynamics of drug resistance. Nyabadza, et al. [9] and White, et al. [10] modeled the dynamics of crystal meth abuse and heroin epidemics. Rwat and co-workers [11] examined the effect of recycling the recovered individuals back into the population while Donoghoe [12] studied the effect of drugs on global health. Murray, et al. (2007) [13] studied the effect of cannabis on mental health Pluddemann, [14] investigated the use of strategies to, monitor alcohol and substance abuse. Recent work of McGinty, et al. [15], Hooker, et al. [16], Butler, et al. [17], Scott, et al. [18] Chang, et al. [19] and Paquette, et al. [20] have demonstrated the devastating effects of drug addiction and the urgent need to combat this problem.

Akanni, et al. [21] Abidemi, et al. [22] and Olaniyi, et al. [23] studied dynamic models involving illicit drug use. All the optimal control work done so far involves single objective minimization. In this work, multiobjective nonlinear model predictive control calculations are performed on drug addiction models in conjunction with bifurcation analysis. It is numerically demonstrated for two problems involving drug addiction that the presence of bifurcation points enables the MNLMPC calculations to converge to the Utopia solution.

The bifurcation analysis and the MNLPMC methods will now be presented followed by an explanation as to why the presence of bifurcation points leads to the MNLMPC calculations converging to the Utopia solution.

Bifurcation analysis

The existence of multiple steady-states (caused by limit and branch point singularities) and oscillatory behavior caused by Hopf bifurcation points) in chemical processes has led to a lot of computational work to explain the causes of these nonlinear phenomena. N MATCONT, [24,25] is a commonly used software to find limit points, branch points, and Hopf bifurcation points. Consider an ODE system.

$\dot{x} = f (x, β) (1)$

$x \in R^{n}$ The tangent plane at any point x is [v₁, v₂, v₃, v₄, …. v_n+1]. Define matrix A given by

$A = [\begin{array}{l} \frac{\partial f_{1}}{\partial x_{1}} \frac{\partial f_{1}}{\partial x_{2}} \frac{\partial f_{1}}{\partial x_{3}} \frac{\partial f_{1}}{\partial x_{4}} .......... \frac{\partial f_{1}}{\partial x_{n}} \frac{\partial f_{1}}{\partial β} \\ \frac{\partial f_{2}}{\partial x_{1}} \frac{\partial f_{2}}{\partial x_{2}} \frac{\partial f_{2}}{\partial x_{3}} \frac{\partial f_{2}}{\partial x_{4}} .......... \frac{\partial f_{2}}{\partial x_{n}} \frac{\partial f_{2}}{\partial β} \\ ........................................................... \\ ........................................................... \\ \frac{\partial f_{n}}{\partial x_{1}} \frac{\partial f_{n}}{\partial x_{2}} \frac{\partial f_{n}}{\partial x_{3}} \frac{\partial f_{n}}{\partial x_{4}} .......... \frac{\partial f_{n}}{\partial x_{n}} \frac{\partial f_{n}}{\partial β} \end{array}] (2)$

With β the bifurcation parameter. The matrix A can be written in a compact form as

$A = [B | \frac{\partial f}{\partial β}] (3)$

The tangent surface must be satisfied.

$A v = 0 (4)$

For both limit and branch points the matrix B must be singular. For a limit point (LP) the n+1th component of the tangent vector v_n+1 = 0 and for a Branch Point (BP) the matrix $[\begin{array}{l} A \\ v^{T} \end{array}]$ must be singular., The function $\det (2 f_{x} (x, β) ⊙ I_{n})$ should be zero for a Hopf bifurcation point. $⊙$ Indicates the alternate product while I_n is the n-square identity matrix. A detailed derivation can be found in Kuznetsov [26,27] and Govaerts [28]. Sridhar [29] used Matcont to perform bifurcation analysis on chemical engineering problems.

MNLMPC (Multiobjective Nonlinear Model Predictive Control) method

The Multiobjective Nonlinear Model Predictive Control (MNLMPC) method was first proposed by Flores Tlacuahuaz, et al. [30] and used by Sridhar [31]. This method is rigorous and it does not involve the use of weighting functions nor does it impose additional parameters or additional constraints on the problem unlike the weighted function or the epsilon correction method [32]. For a problem that is posed as

$\begin{array}{l} \min J (x, u) = (x_{1}, x_{2} .... x_{k}) \\ s u b j e c t t o \frac{d x}{d t} = F (x, u); h (x, u) \leq 0; x^{L} \leq x \leq x^{U}; u^{L} \leq u \leq u^{U} (5) \end{array}$

The MNLMPC method first solves dynamic optimization problems independently minimizing/maximizing each x_i individually. The minimization/maximization of x_i will lead to the values x_i* . Then the optimization problem that will be solved is

$\begin{array}{l} \min \sqrt{{x_{i} - x_{i}^{*}}^{2}} \\ s u b j e c t t o \frac{d x}{d t} = F (x, u); h (x, u) \leq 0; x^{L} \leq x \leq x^{U}; u^{L} \leq u \leq u^{U} (6) \end{array}$

This will provide the control values for various times. The first obtained control value is implemented and the remaining is discarded. This procedure is repeated until the implemented and the first obtained control values are the same.

The optimization package in Python, Pyomo [33] where the differential equations are automatically converted to a Nonlinear Program (NLP) using the orthogonal collocation method [34] is commonly used for these calculations. State-of-the-art solvers like IPOPT [35] and BARON [36] are normally used in conjunction with PYOMO.

Effect of singularities (Limit Point (LP) and Branch Point (BP)) on MNLMPC

Let the minimization be of the variables P₁, P₂ l result in the values M₁ and M₂. The multiobjective objective function to be minimized will result in the problem.

$\begin{array}{l} \min \int_{0}^{t_{f}} {((p_{1} (t) d t) - M_{1})}^{2} + \int_{0}^{t_{f}} {((p_{2} (t) d t) - M_{2})}^{2} = \int_{0}^{t_{f}} Ρ (x, t) d t (7) \\ s u b j e c t t o \frac{d x_{i}}{d t} = g_{i} (x, u) \end{array}$

The Euler Lagrange equation also known as costate equations will be

$\frac{d λ_{i}}{d t} = - (\frac{\partial P}{\partial x_{i}} + λ_{i} g_{i}) (8)$

λ_i is the Lagrangian multiplier. Taking the derivative of the objective function we get

$\frac{d}{d x_{i}} ({(p_{1} - M_{1})}^{2} + {(p_{2} - M_{2})}^{2}) = 2 (p_{1} - M_{1}) \frac{d}{d x_{i}} (p_{1} - M_{1}) + 2 (p_{2} - M_{2}) \frac{d}{d x_{i}} (p_{2} - M_{2}) (9)$

At the Utopia point both (P₁ – M₁) and (P₂ – M₁) are zero. Hence

$\frac{d}{d x_{i}} ({(p_{1} - M_{1})}^{2} + {(p_{2} - M_{2})}^{2}) = 0 (10)$

The co-state equation in optimal control is

$\begin{array}{l} \frac{d}{d t} (λ_{i}) = - \frac{d}{d x_{i}} ((^{p_{1} - M_{1}) 2} + (^{p_{2} - M_{2}) 2}) - g_{x} λ_{i} (11) \\ λ_{i} (t_{f}) = 0 \end{array}$

λ_i is the Lagrangian multiplier. The first term in this equation is 0 and hence.

$\begin{array}{l} \frac{d}{d t} (λ_{i}) = - g_{x} λ_{i} (12) \\ λ_{i} (t_{f}) = 0 \end{array}$

If the set of ODE $\frac{d x}{d t} = g (x, u)$ has a limit or a branch point, g_x is singular.

This implies that there are two different vectors-values for [λ_i] where $\frac{d}{d t} (λ_{i}) > 0$ and $\frac{d}{d t} (λ_{i}) < 0$ . In between there is a vector [λ_i] where This coupled with the boundary condition λ_i(t_f) = 0 will lead to [λ_i] = 0 which will make the problem an unconstrained optimization problem. The only solution for the unconstrained problem is the Utopia solution.

Results and discussion

In this section, the results of bifurcation analysis and MNLMPC calculations for two problems involving drug addiction are presented. The models used are described in Islam, et al. [7] and Mushayabasa, et al. [4]. The equations for each problem are presented followed by the bifurcation analysis and MNLMPC results.

Problem 1 Islam, et al. [7]

Equations representing problem 1

S_a(t) represents individuals who are not drug users, but are at a high risk of taking drugs
L(t)) represents light drug users
H(t) represents heavy drug users
R_v(t) represents drug users under treatment in rehabilitation
Q(t) represents individuals who will never take drugs

The equations are

$\begin{array}{l} \frac{d S_{a}}{d t} = r - α S_{a} H - μ S_{a} - u_{1} S_{a} \\ \frac{d L}{d t} = α S_{a} H - μ L - β L - δ L - u_{2} L \\ \frac{d H}{d t} = β L - μ H - γ H + p_{a} R - u_{3} H (13) \\ \frac{d R v}{d t} = γ H - μ R - θ R - p_{a} R \\ \frac{d Q}{d t} = θ R - μ Q + u_{1} S_{a} + u_{2} L + u_{3} H \end{array}$

The model parameters are

r = 4.25; µ = 0.00561; α = 0.002; β = 0.6; δ = 0.025; γ = 1.5; pa = 0.02

u₁,u₂,u₃ are the control variables

where

r represents the recruitment rate of the population
µ is the natural mortality rate
α is the interaction rate among the susceptible and light drug users
β is the effective rate at which light users convert into heavy drug users
δ the removal rate from addiction without treatment
γ is the rate at which heavy addicts are being sent to rehabilitation for treatment
u₁ is the awareness and educational programs
u₂ is the family-based care
u₃ represents the effectiveness of rehabilitation centers

Bifurcation analysis for problem 1

When bifurcation analysis with µ the bifurcation parameter was performed on the equations representing problem 1, a branch point was found a [S_a, L, H, R_v, Q, µ] values of (782.26, 0.0, 0, 0,0, 0.005433 ). Figure 1a shows the bifurcation diagram with this branch point.

MLNMPC for problem 1

The MNLMPC of problem 1, $\sum Q (t)$ was maximized and resulted in a value of 2000; while $\sum H (t)$ was minimized and resulted in a value of 0. The multiobjective optimal control problem involved the minimization of $(\sum Q (t) - 2000)^{2} + (\sum H (t) - 0)^{2}$ the subject to the dynamic equation set representing this problem. This resulted in the Utopia point of 0 and the MNLMPC values of the control variables obtained were [u₁, u₂, u₃] = [0.0004, 0.0405, 0.5362]. The MNLMPC profiles are shown in Figures 1a-1i.

Problem 2 Mushayabasa, et al. [4]

Equations representing problem 2

In this problem, the time-dependent variables are

S_v(t) susceptible individuals
I(t) light or occasional drug users
I_av(t) heavy drug users
M_v(t) mentally ill population and (individuals who suffer mental illness due to drug use,
R_v(t) detected illicit drug users

The equations that represent the drug addiction problem are

$\begin{array}{l} \frac{d S_{v}}{d t} = μ - (1 - u_{c}) λ S_{v} - μ S_{v} \\ \frac{d I_{v}}{d t} = (1 - u_{c}) λ S_{v} - (α + γ v_{c} + σ + μ + ψ) I_{v} (14) \\ \frac{d I_{a v}}{d t} = α I_{v} - (ρ v_{c} + ϕ + μ + d) I_{a v} \\ \frac{d M_{v}}{d t} = σ I_{v} + ϕ I_{a v} - (μ + ε v_{c} + δ) M_{v} \\ \frac{d R_{v}}{d t} = v_{c} (γ I + ρ I_{a v} + ε M_{v}) - (μ + ω) R_{v} \\ λ = β (I_{v} + k I_{a v}) \end{array}$

And the parameter values are,

ω = 0.3; µ = 0.02; k = 1.25; β = 0.35; γ = 0.1; ρ = 0.35; ε = 0.6; α = 0.01; ψ = 0.03; δ = 0.14; d = 0.2; σ= 0.05;ϕ = 0.09;

uc, vc are the control variables

Here,

α represents the rate at which light drug users become heavy drug users
γ, ε, ρ the rates of detection and rehabilitation of individuals in classes I_v, M_v, I_av
σ, φ the rates at which light and heavy illicit drug users develop mental illness
ψ, d the permanent exit rates of light and heavy users
δ mentally ill individuals who permanently exit the model because of death
ω the rate at which individuals recover as a result of rehabilitation
β the strength of interactions between susceptible individuals and illicit drug users
u_c represents the reduction of the intensity of “social influence”
v_c models the effort on the detection of illicit drug users

Bifurcation analysis for problem 2

When bifurcation analysis with $α$ as the bifurcation parameter was performed on the equations representing problem 2, a branch point was found at [S_v, I_v, I_av, M_v, R_v, α] = [1.0, 0.0, 0.0, 0.0, 0.0, 0.430112]. The bifurcation diagram is shown in Figure 2a.

MLNMPC for problem 2

For the MNLMPC of problem 2, $\sum I_{v} (t)$ and $\sum I_{a v} (t)$ were minimized individually and both the minimizations resulted in a value of 0. The multiobjective optimal control problem involved the minimization of $(\sum I_{v} (t))^{2} + (\sum I_{a v} (t))^{2}$ the subject to the dynamic equation set representing this problem. This resulted in the Utopia point of 0 and the MNLMPC values of the control variables obtained were [u₁, u₂, u₃] = [0.0004, 0.0405, 0.5362]. The various MNLMPC profiles are shown in Figures 2b-2h.

Two problems involving drug addiction models have been shown to exhibit branch points leading to two different solution branches. In both cases, it is computationally shown that the MNLMPC calculations would converge to the Utopia solution as the theoretical analysis predicts.

Conclusions and future work

Branch points leading to two separate branches were exhibited when bifurcation analysis was performed on the two drug addiction models considered in this paper. The rigorous analysis demonstrated that the presence of the branch points would result in the MNLMPC calculations converging to the Utopia solution. This fact was also computationally validated. The presence of the branch points indicates that there are two paths that the drug abuse dynamics can take. The attainment of the Utopia point shows that the number of light and heavy drug users can both be minimized while maximizing the number of users who do not take drugs. Future work would involve using drug addiction models with time delay.