Marc Bigirimana^1*, Fulgence Nahayo¹, Mounir Haddou², Abraham Niyongere^3
1Laboratoire Universitaire de Recherche en Modélisation et en Ingénierie de Statistique Appliquée, Université du Burundi, Bujumbura, Burundi.
²Institut de Recherche Mathématique de Rennes/Institut National des Sciences Appliquées de Rennes, Université de Rennes, Rennes, France.
³Centre de Recherche en Infrastructures, Environnement et Technologie, Université du Burundi, Bujumbura, Burundi.
DOI: 10.4236/ojapps.2026.167134   PDF    HTML   XML   4 Downloads   42 Views   Citations

Abstract

The application of Pontryagin’s Maximum Principle (PMP) to optimize the dynamics of vehicle queues in road traffic systems. We present a mathematical framework that integrates optimal control theory with traffic flow dynamics to enhance the efficiency and performance of road networks. By modeling vehicle behavior and interactions within queues, we identify optimal control strategies that minimize delays and maximize throughput. The proposed methodology includes a comprehensive analysis of system constraints and environmental factors impacting traffic flow. Simulation results demonstrate the effectiveness of the PMP approach in various traffic scenarios, highlighting significant improvements in queue management and overall traffic efficiency. Our findings contribute to the development of advanced traffic management systems and provide valuable insights for future research in transportation optimization.

Keywords

State Constraint Penalization, Optimal Control Problem, Exterior Penalization, Traffic Fluidity Function, Traffic Queue Optimization, Joint Optimization Problem, Pontryagin Maximum Principle, Symplectic Discretization

Share and Cite:

1. Introduction

Efficient management of road traffic is a central challenge in modern transportation systems, with direct implications for travel time, fuel consumption, and environmental impact. Vehicle queues at bottlenecks intersections, on-ramps, signalized junctions, or lane drops are a primary source of delay and flow degradation. This paper investigates a control-theoretic approach to reducing queueing delay and improving throughput by applying Pontryagin’s Maximum Principle (PMP) to a penalized joint optimization problem for vehicle-queue dynamics.

In particular, we show that a coupled two-vehicle continuous-time dynamics model can serve as a building block for representing queue evolution: the interaction between the lead and the following vehicles captures the local propagation of congestion, while the combined state and control trajectories approximate the queue’s aggregate behavior.

We construct a continuous-time optimal control model that captures individual vehicle interactions within a platoon and the aggregate behavior of the queue, incorporate operational and safety constraints, and introduce penalty terms to enforce feasibility and desirable performance metrics. Using the indirect (PMP-based) approach, we derive necessary conditions for optimality, obtain a reduced Hamiltonian formulation where possible, and develop a symplectic Partitioned Runge-Kutta (SPRK) discretization to integrate the state-adjoint system while preserving geometric structure. Numerical experiments and simulations illustrate how the PMP framework yields control laws that reduce delays and smooth traffic flow under various scenarios. The contributions of this work are: i) a rigorous PMP formulation for joint optimization of vehicle queues with penalty terms for constraints; ii) an analysis of regularity and implementability conditions enabling the implicit elimination of control variables; and iii) a structure-preserving numerical scheme (SPRK) adapted to the resulting Hamiltonian boundary-value problem, validated with simulation results.

2. Literature Review

2.1. Optimal Control and Traffic Flow

The intersection of optimal control theory and traffic modeling has produced a rich literature addressing ramp metering, signal timing, and ramp-to-mainline coordination (See Papageorgiou et al., 2008 [1]; Daganzo, 1994 [2]). Classical macroscopic models (Lighthill-Whitham-Richards, kinematic waves) and microscopic car-following models have both been used as the dynamical substrate for control design. More recent works concentrate on platoon control and cooperative adaptive cruise control (CACC) where local interactions are optimized to improve string stability and throughput (Swaroop & Hedrick, 1996 [3]; Rajamani, 2012 [4]).

2.2. Pontryagin’s Maximum Principle in Traffic Applications

PMP has been applied to various transportation problems including speed profile optimization, eco-driving, and gap-minimizing strategies for vehicle merging (See Miculescu & Karaman, 2015 [5]; Zhao, L., Malikopoulos et al., 2018 [6]). These works exploit PMP’s capability to provide necessary optimality conditions and to characterize bang-bang or singular controls in constrained settings. The indirect method based on PMP is particularly valuable when analytical insight on the control structure is sought, although it requires solving two-point boundary value problems (TPBVPs) and ensuring appropriate regularity.

2.3. Penalization and Constrained Optimal Control

Handling state and control constraints in PMP formulations often involves penalization or the use of Lagrange multipliers and complementarity conditions. Classical texts (See, Bryson & Ho, 2018 [7]; White, David H., 2012 [8]) detail penalization techniques and constraint relaxation methods. More recent computational approaches combine penalty terms with interior point solvers or augmented Lagrangian methods to tractably enforce constraints in large-scale control problems (See Nocedal & Wright, 2006 [9]; Wächter & Biegler, 2006 [10]).

2.4. Discrete Structure-Preserving Integrators for Hamiltonian Systems

When using the indirect PMP approach the resulting necessary conditions form a Hamiltonian boundary-value problem. Symplectic integrators, including symplectic Runge-Kutta and Partitioned Runge-Kutta (PRK) methods, have been shown to preserve qualitative features of Hamiltonian flows (hairiness: conservation of the symplectic form, near-conservation of energy over long times) and to yield improved long-term numerical behavior for TPBVPs (Hairer, Lubich & Wanner, 2006 [11]; Sanz-Serna & Calvo, 1995 [12]). The theory of B-series and P-series provides a systematic way to derive order conditions for such methods (Butcher, 1987 [13]; Laurent-Varin, J. 2005 [14]).

2.5. Applications of PMP with Structure-Preserving Discretization in Traffic

A growing body of work addresses the numerical implementation of PMP for traffic problems with attention to preserving system structure. Examples include optimal control of connected automated vehicle platoons using indirect methods with tailored integrators (Hairer, E., Lubich et al., 2002 [15]) and trajectory optimization for eco-driving with symplectic or variational integrators to maintain stability over long horizons (Zhang, Y., & Cassandras, C. G. 2019 [16]).

2.6. Gaps and Motivation for the Present Study

Despite extensive literature on PMP and traffic control, relatively few studies explicitly combine a penalized joint optimization for vehicle queues with a rigorous symplectic discretization tailored to the resulting Hamiltonian structure. In particular, the treatment of multi-vehicle coupling, inequality constraints on states and controls, and the numerical preservation of geometric invariants in the discrete TPBVP setting remains underexplored. This paper aims to fill that gap by: formulating a penalized PMP model for joint queue optimization, analyzing regularity conditions that permit control elimination, and devising a practical SPRK4 discretization whose coefficients satisfy the required order and symplecticity conditions, validated through AMPL/IPOPT and MATLAB simulations.

3. Mathematical Analysis of Penalized Joint Optimization of Vehicle Queue Dynamics in Road Traffic

Consider now an optimal control problem with constraints of the following form:

$\{\begin{array}{l}\underset{\delta \in \mathcal{U}}{\text{min}}J\left(y(.),\delta (.)\right)\hfill \\ \dot{y}\left(t\right)=f\left(y\left(t\right),\delta \left(t\right)\right)\hfill \\ g^L\le g\left(y,\delta \right)\le g^U\hfill \\ y\left(t_0\right)=y_0\hfill \end{array}$(1)

In system (1), $f:ℝ\times {ℝ}^n\times {ℝ}^m\to {ℝ}^n$ and $g:{ℝ}^n\times {ℝ}^m\to {ℝ}^n$ are maps of class $C^{\infty }$ . The functions $g^L$ and $g^U$ denote, respectively, the lower and upper bounds on the state variables.

Penalty Techniques for Constraints in Optimal Control

In the context of optimal control problems, it is often necessary to address constraints on states or controls, which can complicate the problem’s resolution. Penalty methods provide an effective approach to integrate these constraints into the optimization process. This section focuses on two types of state constraint penalization, exterior penalization and interior penalization. Exterior penalization involves adding a penalty term to the objective function to discourage solutions that violate the constraints. This method is particularly useful when the constraints are difficult to manage or when adherence to them is critical [17]. So, the penalized objective function in the system (1) becomes:

$J_{\text{pen}}(.)=J(.)+\frac{\rho }{2}{\displaystyle \sum }_{i=1}^n\left\{{\left(\text{max}\left(0,g_i\left(y,\delta \right)-g_i^U\right)\right)}^2+{\left(\text{max}\left(0,g_i^L-g_i\left(y,\delta \right)\right)\right)}^2\right\},$ (2)

where $\rho >0$ is a penalty parameter that controls the significance of the constraint. It should be noted that in a problem of the form (1), the exterior penalty affects only the objective function $J(.)$ and not the constraints $g^L\le g\left(y,\delta \right)\le g^U$ or the initial conditions $y\left(t_0\right)=y_0$ .

On the other hand, interior or logarithmic penalization aims to keep the solutions within the constraint space by modifying the cost function. This approach can be more effective when the constraints are defined to allow admissible solutions at the boundaries.

$J_{\text{pen}}(.)=J(.)-\epsilon {\displaystyle \sum }_{i=1}^n\left\{\text{ln}\left(g_i^U-g_i\left(y,\delta \right)\right)+\text{ln}\left(g_i\left(y,\delta \right)-g_i^L\right)\right\}$ (3)

The logarithmic term in (3) penalizes solutions that approach the boundaries of the constraints, and the scalar parameter $\epsilon$ must be small. In this specific case, an admissible solution to this problem is such that the state equation $\dot{y}=f\left(y,\delta \right)$ remains valid, as well as the inequalities $g_i^U-g_i\left(y,\delta \right)>0$ and $g_i\left(y,\delta \right)-g_i^L>0$ [14] [18]. Penalty methods, whether exterior or interior, offer flexible approaches to incorporate state constraints into optimal control problems. The choice between these methods often depends on the nature of the constraints, the characteristics of the problem, and the desired performance criteria. However, selecting the penalty parameter is an iterative process that may require experimentation.

4. SPRK Discretization of a Constrained Penalized Optimal Control Problem

4.1. General Description of Symplectic PRK Methods

The general principle of SPRK methods is based on the system dynamics written as a system of differential equations of the form:

$\dot{y}=f\left(y,\delta \right),y\left(0\right)=y_0,$(4)

which can be rewritten in the form of the following Hamiltonian system

$\dot{y}={\tilde{J}}^{-1}\nabla \mathscr{H}\left(y,\tilde{p}\right),$ (5)

where $\tilde{J}=\left(\begin{array}{cc}0& I_n\\ I_n& 0\end{array}\right)$ is the canonical matrix for Hamiltonian systems and $I_n$ denotes the $n\times n$ identity matrix [19]. Note that, for the Hamiltonian system (5), the following Jacobian matrix of partial derivatives of the numerical flow with respect to $y$ and $\tilde{p}$ defined by:

$\begin{array}{lll}{\xi }_h:{ℝ}^{2n}\hfill & \to \hfill & {ℝ}^{2n}\hfill \\ \left(y_n,{\tilde{p}}_n\right)\hfill & \mapsto \hfill & \left(y_{n+1},{\tilde{p}}_{n+1}\right)\hfill \end{array}$

preserves the canonical symplectic form and thus satisfies the (6) following Poincaré matrix identity [20]

${\left(\frac{\partial {\xi }_h}{\partial \left(y,\tilde{p}\right)}\right)}^{\text{T}}\tilde{J}\left(\frac{\partial {\xi }_h}{\partial \left(y,\tilde{p}\right)}\right)=\tilde{J},$(6)

then the Runge-Kutta numerical method $y_{n+1}={\xi }_h\left(y_n\right)$ is said to be symplectic for the Hamiltonian (5), with $\mathscr{H}$ regular and for every sufficiently small time step [19] [21]. We say that a PRK scheme is symplectic if its associated flow is symplectic. It is known that PRK schemes satisfying the (7) following relation:

${\widehat{b}}_i:=b_i,{\widehat{a}}_{ik}:=b_k-\frac{b_k}{b_i}{a}_{ki},i=1,\cdots ,s\text{;}k=1,\cdots ,s.$(7)

are symplectic; we denote this class by SPRK [15] [22] [23].

4.2. Optimal Control Problem Formulation

The complete statement of the optimal control problem addressed in this work, together with its detailed description, is presented using the indirect approach based on the Pontryagin Maximum Principle. This approach stems from the joint optimization of a platoon of road vehicles, as it is described in [24]. We implement and analyze this method here to derive the necessary optimality conditions and to develop an appropriate symplectic discretization.

Before extending the framework to a full queue of vehicles, the interaction dynamics are first studied at the pairwise level, considering a leader-follower pair as the elementary unit of analysis. This two-vehicle setting allows for a rigorous examination of the local coupling between consecutive vehicles, providing the foundational understanding necessary for scaling the model to larger platoons.

In the following, the complete problem formulation is:

$\{\begin{array}{l}\underset{\left(y,\delta \right)\in {ℝ}_{ad}^{15}\times {ℝ}_{ad}^{15}\times \mathcal{U}}{\text{max}}{J}_{G_{12}}\left(y(.),\delta (.)\right)=\lambda J\left(y_1\left(t\right),{\delta }_1\left(t\right)\right)+\left(1-\lambda \right)J\left(y_2\left(t\right),{\delta }_2\left(t\right)\right)+P(.)\\ \dot{y}\left(t\right)=f\left(y\left(t\right),\delta \left(t\right)\right),y\left(t\right)=\left(y_1\left(t\right),y_2\left(t\right)\right),\delta \left(t\right)=\left({\delta }_1\left(t\right),{\delta }_2\left(t\right)\right)\\ l_{1j}\left(y_j\left(t\right),{\delta }_j\left(t\right)\right)\le 0,l_{2j}\left(y_j\left(t\right),{\delta }_j\left(t\right)\right)\ge 0\\ y_{j_{\text{min}}}\le y_j\left(t\right)\le y_{j_{\text{max}}},{\delta }_{j_{\text{min}}}\le {\delta }_j\left(t\right)\le {\delta }_{j_{\text{max}}},j=1,2\end{array}$ (8)

In the formulation (8) above, $J_{G_{12}}(.)$ is the coupled-model cost; i.e. the overall smoothness, $\dot{y}$ denotes the traffic dynamics for a pair of vehicles, while $l_{1j}\left(y_j,{\delta }_j\right)$ and $l_{2j}\left(y_j,{\delta }_j\right)$ are constraints on the state and the control, respectively. The function $J\left(y_j,{\delta }_j\right)=\frac{1}{v_{\max }^2}{\Vert R\left({\varphi }_j,{\theta }_j,{\psi }_j\right)v_{a_j}+{\Omega }_j\times r_{G_j}^o\Vert }^2+g\left(y_j,{\delta }_j\right)$ is the individual traffic smoothness function, where $R\left({\varphi }_j,{\theta }_j,{\psi }_j\right)$ is the rotation operator matrix, $v_{a_j}={\left(\begin{array}{l}{u}_j,v_j,w_j\hfill \end{array}\right)}^{\text{T}}$ is the individual vehicle aerodynamic velocity, ${\Omega }_j={\left(\begin{array}{l}{p}_j,q_j,r_j\hfill \end{array}\right)}^{\text{T}}$ is the vehicle’s angular velocity, $r_{G_j}^o={\left(\begin{array}{l}{X}_{G_j}^o,Y_{G_j}^o,Z_{G_j}^o\hfill \end{array}\right)}^{\text{T}}$ is the vehicle’s vector position in the $R_{G_j}$ frame relative to the observer’s $R_O$ frame, $y_j$ is the vehicle’s state vector and $g\left(y_j,{\delta }_j\right)=\frac{{\Vert {\delta }_j\Vert }^2}{{\delta }_{j\text{max}}^2}$ is the control price function. The term $P(.)={\eta }_1\left[\frac{{\left(S_{G_{12}}^o-S_{\text{min}}\right)}^2}{1+{\left(S_{G_{12}}^o-S_{\text{min}}\right)}^2}\right]+{\eta }_2\left[\frac{{\left({\left(Y_{G_{12}}^{o^2}+{\epsilon }^2\right)}^{1/2}-Y_{\text{max}}\right)}^2}{1+{\left({\left(Y_{G_{12}}^{o^2}+{\epsilon }^2\right)}^{1/2}-Y_{\text{max}}\right)}^2}\right]+{\eta }_3\left(u_{12}^2+v_{12}^2\right)$ is the penalty function that includes penalty terms, encouraging minimal separation and appropriate relative speeds [25] [26].

However, in the previous formulation, certain aircraft-inspired variables such as aerodynamic velocity, roll, pitch, yaw, and observer-frame coordinates are relevant because they provide a fine-grained dynamic description of how vehicle states as velocities and orientations evolve over time, enabling a direct link between congestion and the vehicles’physical behavior. The goal is to exploit the propagation of congestion and inter-vehicle perturbations to model the queue’s overall evolution. Accordingly, transferring the dynamics to the queue-level description is justified by approximating the aggregate behavior from leader-follower interactions under operational constraints.

4.3. Optimality Conditions and Symplectic Discretization

For any absolutely continuous function $\tilde{p}(.):\left[0,T\right]\to {ℝ}_0^n$ , called the adjoint state vector, and any real number ${\tilde{p}}_0\le 0$ , such that the pair $\left(\tilde{p}(.),{\tilde{p}}_0\right)$ is nontrivial, the pseudo-Hamiltonian of problem (8) is defined by the following function:

$H\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)=\langle {\tilde{p}}^{\text{T}},f\left(y,\delta ,t\right)\rangle -{\tilde{p}}_0{J}_{G_{12}}\left(y,\delta ,t\right)-{\mu }_1{l}_{1j}\left(y,\delta ,t\right)+{\mu }_2{l}_{2j}\left(y,\delta ,t\right)$ (9)

where $\tilde{p}=\left({\tilde{p}}_1,{\tilde{p}}_1\right)$ and ${\tilde{p}}_j\in {ℝ}^{15}$ is the adjoint vector. With the new Pontryagin formulation (9), we state the necessary optimality conditions as follows:

$\begin{array}{l}{H}_{\delta }\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)=0\\ \dot{y}=H_{\tilde{p}}\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)=f\left(y,\delta ,t\right)\\ \dot{\tilde{p}}=-H_y\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)\\ l_{1j}\left(y,\delta ,t\right)\le 0,l_{2j}\left(y,\delta ,t\right)\ge 0\\ t\in \left[0,T\right],{\mu }_1\le 0,{\mu }_2\ge 0\end{array}$(10)

The state-adjoint equations are derived from:

$\begin{array}{l}{\dot{\tilde{p}}}_{1j}=-\frac{\partial H}{\partial m_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{2j}=-\frac{\partial H}{\partial v_{a_j}},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{3j}=-\frac{\partial H}{\partial {\gamma }_{a_j}},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{4j}=-\frac{\partial H}{\partial u_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{5j}=-\frac{\partial H}{\partial v_j},\\ {\dot{\tilde{p}}}_{6j}=-\frac{\partial H}{\partial w_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{7j}=-\frac{\partial H}{\partial p_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{8j}=-\frac{\partial H}{\partial q_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{9j}=-\frac{\partial H}{\partial r_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{10j}=-\frac{\partial H}{\partial {\theta }_j},\\ {\dot{\tilde{p}}}_{11j}=-\frac{\partial H}{\partial {\varphi }_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{12j}=-\frac{\partial H}{\partial {\psi }_j},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{13j}=-\frac{\partial H}{\partial X_{G_j}^o},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{14j}=-\frac{\partial H}{\partial Y_{G_j}^o},\mathrm{<mi>\; </mi>}{\dot{\tilde{p}}}_{13j}=-\frac{\partial H}{\partial Z_{G_j}^o}\end{array}$ (11)

From the constraint relation $H_{\delta }\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)=0$ , we deduce certain regularity assumptions. Let $\left(\overline{u},\overline{y},\overline{p}\right)$ be an extremal that satisfies the necessary conditions optimality as it’s described in (10). If the map

$\delta \mapsto H_{\delta \delta }\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)$ is invertible along the trajectory,(12)

then by the implicit function theorem, in a small $L^{\infty }$ -neighborhood of that trajectory we have $H_{\delta }\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right)=0$ iff $\delta (.)$ can be expressed as a smooth function of the form $\delta \left(t\right)=\varphi \left(y\left(t\right),p\left(t\right)\right)$ . Under these conditions, the true Hamiltonian of the system is defined as the following:

$\mathscr{H}\left(y,\delta ,\tilde{p},{\tilde{p}}_0,t\right):=H\left(\varphi \left(y,\tilde{p}\right),y,\tilde{p}\right)$ (13)

By using $H_{\delta }\left(\varphi \left(y\left(t\right),\tilde{p}\left(t\right)\right),y\left(t\right),\tilde{p}\left(t\right)\right)=0$ in (13), we have ${\mathscr{H}}_y\left(y,\tilde{p}\right)=H_y\left(\varphi \left(y,\tilde{p}\right),y,\tilde{p}\right)$ and ${\mathscr{H}}_{\tilde{p}}\left(y,\tilde{p}\right)=H_{\tilde{p}}\left(\varphi \left(y,\tilde{p}\right),y,\tilde{p}\right)$ . Therefore, under assumption (12), the optimality conditions in (10) are locally equivalent to the following reduced Hamiltonian system [27]:

$\begin{array}{l}\dot{y}={\mathscr{H}}_{\tilde{p}}\left(y,\tilde{p}\right),y\left(0\right)=y_0\\ \dot{\tilde{p}}=-{\mathscr{H}}_y\left(y,\tilde{p}\right),\tilde{p}\left(T\right)={\varphi }^{\prime }\left(y\left(T\right)\right)\\ l_{1j}\left(y,\delta ,t\right)\le 0,l_{2j}\left(y,\delta ,t\right)\ge 0,t\in \left[0,T\right]\end{array}$(14)

According to Lagrange, a curve $\left(y_n,{\tilde{p}}_n,{\delta }_{n_i},{\tilde{q}}_{n_i}\right)$ is an optimal solution of the (15) following discretized problem:

$\{\begin{array}{l}\underset{\left(y_n,{\delta }_{n_i}\right)}{\max }{J}_{G_{12}}\left(y_n,{\delta }_{n_i}\right)\hfill \\ 0=y_n-y_{n+1}+h{\displaystyle \sum _{i=1}^s}{b}_{i}f\left(y_{n_i},{\delta }_{n_i}\right)\hfill \\ 0=y_n-y_{n_i}+h{\displaystyle \sum _{k=1}^s}{a}_{ik}f\left(y_{n_k},{\delta }_{n_k}\right)\hfill \\ 0=y_0-y\left(0\right)\hfill \end{array}$ (15)

when the five following conditions in (16) are simultaneously satisfied:

$\begin{array}{l}{L}_{y_n}={\nabla }_{y_n}\left(L_n+L_{n+1}\right)=0\\ L_{y_{n_i}}={\nabla }_{y_{n_i}}{L}_{n+1}=0,\mathrm{<mi>\; </mi>}i=1,\cdots ,s\\ L_{{\delta }_{n_i}}={\nabla }_{{\delta }_{n_i}}{L}_{n+1}=0,\mathrm{<mi>\; </mi>}i=1,\cdots ,s\\ L_{y_{{\tilde{p}}_n}}={\nabla }_{{\tilde{p}}_n}{L}_n=0,\mathrm{<mi>\; </mi>}n=0,\cdots ,N_t-1\\ L_{y_{{\tilde{q}}_{n_i}}}={\nabla }_{{\tilde{q}}_{n_i}}{L}_{n+1}=0,\mathrm{<mi>\; </mi>}i=1,\cdots ,s\end{array}$ (16)

where $h$ is the discretization step size, $A=\left(a_{ik}\right)$ , and $b=\left(b_i\right)$ are Runge-Kutta coefficients. The Lagrangian function associated with (15) reads:

$\begin{array}{c}L=J_{G_{12}}\left(y_n,{\delta }_{n_i}\right)+{\tilde{p}}_0^{\text{T}}\left(y_0-y\left(0\right)\right)+{\displaystyle \sum }_{n=0}^{N_t-1}\{{\tilde{p}}_{n+1}^{\text{T}}\cdot \left(y_n-y_{n+1}+h{\displaystyle \sum }_{i=1}^s{b}_{i}f\left(y_{n_i},{\delta }_{n_i}\right)\right)\\ +{\displaystyle \sum }_{i=1}^s{\tilde{q}}_{n_i}^T\cdot \left(y_n-y_{n_i}+h{\displaystyle \sum }_{k=1}^s{a}_{ik}f\left(y_{n_k},{\delta }_{n_k}\right)\right)\}\end{array}$ (17)

${\tilde{p}}_{n+1}$ , ${\tilde{q}}_{n_i}$ and ${\tilde{p}}_0$ are Lagrange multipliers associated with system (15). In what follows, the variables, ${\tilde{p}}_n$ will be interpreted as the discretization of the adjoint state from the continuous formulation. We have the following optimality conditions:

$\begin{array}{l}{\tilde{p}}_{N_t}={\Phi }^{\prime }\left(y_{N_t}\right),\mathrm{<mi>\; </mi>}{\tilde{p}}_1={\tilde{p}}_0,\\ {\tilde{p}}_n-{\tilde{p}}_{n+1}={\displaystyle \sum _{i=1}^s}{f}_y{\left(y_{n_i},{\delta }_{n_i}\right)}^{\text{T}}{\tilde{q}}_{n_i},\\ 0=h{b}_i{\tilde{p}}_{n+1}+h{\displaystyle \sum _{k=1}^s}{a}_{ki}{f}_y{\left(y_{n_k},{\delta }_{n_k}\right)}^{\text{T}}{\tilde{q}}_{n_k}-{\tilde{q}}_{n_i},\\ 0=f_{\delta }{\left(y_{n_i},{\delta }_{n_k}\right)}^{\text{T}}{\tilde{q}}_{n_i},\mathrm{<mi>\; </mi>}n=0,\cdots ,N_t-1;\mathrm{<mi>\; </mi>}i=1,\cdots ,s.\end{array}$(18)

Using the assumption that $b_i\ne 0$ and defining ${\tilde{p}}_{n_i}:={\tilde{q}}_{n_i}/\left(h{b}_i\right)$ , $\forall n=0,\cdots ,$ $N_t-1$ and $i=1,\cdots ,s$ , the elimination of ${\tilde{q}}_{n_i}$ in (18) leads to the application of a Partitioned Symplectic Runge-Kutta method of types $M=\left(A,b\right)$ and $\widehat{M}=\left(\widehat{A},\widehat{b}\right)$ , which yields the discrete optimality conditions of the following form:

$\{\begin{array}{l}{y}_{n+1}=y_n+h{\displaystyle \sum _{i=1}^s}{b}_{i}f\left(y_{n_i},{\delta }_{n_i}\right),\mathrm{<mi>\; </mi>}y\left(0\right)=y_0\hfill \\ y_{n_i}=y_n+h{\displaystyle \sum _{k=1}^s}{a}_{ik}f\left(y_{n_k},{\delta }_{n_k}\right)\hfill \\ {\tilde{p}}_{n+1}={\tilde{p}}_n-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_y\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{N_t}={\Phi }^{\prime }\left(y_{N_t}\right)\hfill \\ {\tilde{p}}_{n_i}={\tilde{p}}_n-h{\displaystyle \sum _{k=1}^s}{\widehat{a}}_{ik}{H}_y\left(y_{n_k},{\delta }_{n_k},{\tilde{p}}_{n_k}\right)\hfill \\ 0=H_{\delta }\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)\hfill \end{array}$(19)

4.4. Numerical Implementation of the SPRK4 Algorithm

The SPRK4 algorithm proceeds as follows:

1) We discretize the time interval $\left[0,T\right]$ into $N_t$ subintervals of step size $h=t_{n+1}-t_n=\frac{T}{N_t}$ , where $N_t$ is the maximum number of iterations.

2) For $0\le n\le N_t$ :

$\begin{array}{l}{y}_{n+1}=y_n+h{\displaystyle \sum _{i=1}^s}{b}_{i}f\left(y_{n_i},{\delta }_{n_i}\right),\mathrm{<mi>\; </mi>}y\left(0\right)=y_0\\ y_{n_i}=y_n+h{\displaystyle \sum _{k=1}^s}{a}_{ik}f\left(y_{n_k},{\delta }_{n_k}\right)\\ {\tilde{p}}_{n+1}={\tilde{p}}_n-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_y\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{N_t}={\Phi }^{\prime }\left(y_{N_t}\right)\\ {\tilde{p}}_{n_i}={\tilde{p}}_n-h{\displaystyle \sum _{k=1}^s}{\widehat{a}}_{ik}{H}_y\left(y_{n_k},{\delta }_{n_k},{\tilde{p}}_{n_k}\right)\\ H_{\delta }\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)=0\\ {\widehat{a}}_{ik}:=b_k-\frac{b_k}{b_i}{a}_{ki},\mathrm{<mi>\; </mi><mi>\; </mi>}{\widehat{b}}_i=b_i\end{array}$(20)

3) Write: $t_{n+1}$ , $y_{n+1}$ , ${\tilde{p}}_{n+1}$ .

4) Stop

The first equation of (20) can be written in the following explicit form:

(21)

The second equation of (20) can be written in the following explicit form:

(22)

The third equation of (20) can be written in the following explicit form:

$\begin{array}{l}{\tilde{p}}_{1j,n+1}={\tilde{p}}_{1j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{m_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{1j,N_t}={\tilde{p}}_{1j}\left(T\right)\\ {\tilde{p}}_{2j,n+1}={\tilde{p}}_{2j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{v_{aj}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{2j,N_t}={\tilde{p}}_{2j}\left(T\right)\\ {\tilde{p}}_{3j,n+1}={\tilde{p}}_{3j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{{\gamma }_{aj}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{3j,N_t}={\tilde{p}}_{3j}\left(T\right)\\ {\tilde{p}}_{4j,n+1}={\tilde{p}}_{4j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{u_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{4j,N_t}={\tilde{p}}_{4j}\left(T\right)\\ {\tilde{p}}_{5j,n+1}={\tilde{p}}_{5j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{v_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{5j,N_t}={\tilde{p}}_{5j}\left(T\right)\\ {\tilde{p}}_{6j,n_i}={\tilde{p}}_{6j,n}-h{\displaystyle \sum _{k=1}^s}{\widehat{a}}_{ik}{H}_{w_j}\left(y_{n_k},{\delta }_{n_k},{\tilde{p}}_{n_k}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}i=1,\cdots ,s\end{array}$

$\begin{array}{l}{\tilde{p}}_{7j,n+1}={\tilde{p}}_{7j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{p_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{7j,N_t}={\tilde{p}}_{7j}\left(T\right)\\ {\tilde{p}}_{8j,n+1}={\tilde{p}}_{8j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{q_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{8j,N_t}={\tilde{p}}_{8j}\left(T\right)\\ {\tilde{p}}_{9j,n+1}={\tilde{p}}_{9j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{r_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{9j,N_t}={\tilde{p}}_{9j}\left(T\right)\\ {\tilde{p}}_{10j,n+1}={\tilde{p}}_{10j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{{\theta }_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{10j,N_t}={\tilde{p}}_{10j}\left(T\right)\\ {\tilde{p}}_{11j,n+1}={\tilde{p}}_{11j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{{\varphi }_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{11j,N_t}={\tilde{p}}_{11j}\left(T\right)\\ {\tilde{p}}_{12j,n+1}={\tilde{p}}_{12j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{{\psi }_j}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{12j,N_t}={\tilde{p}}_{12j}\left(T\right)\\ {\tilde{p}}_{13j,n+1}={\tilde{p}}_{13j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{X_{G_j}^o}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{13j,N_t}={\tilde{p}}_{13j}\left(T\right)\\ {\tilde{p}}_{14j,n+1}={\tilde{p}}_{14j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{Y_{G_j}^o}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{14j,N_t}={\tilde{p}}_{14j}\left(T\right)\\ {\tilde{p}}_{15j,n+1}={\tilde{p}}_{15j,n}-h{\displaystyle \sum _{i=1}^s}{\widehat{b}}_i{H}_{Z_{G_j}^o}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right),{\tilde{p}}_{15j,N_t}={\tilde{p}}_{15j}\left(T\right)\end{array}$(23)

The fourth equation of (20) can be written in the following explicit form:

$\begin{array}{l}{\tilde{p}}_{14j,n_i}={\tilde{p}}_{14j,n}-h{\displaystyle \sum _{k=1}^s}{\widehat{a}}_{ik}{H}_{Y_{G_j}^o}\left(y_{n_k},{\delta }_{n_k},{\tilde{p}}_{n_k}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}i=1,\cdots ,s\\ {\tilde{p}}_{15j,n_i}={\tilde{p}}_{15j,n}-h{\displaystyle \sum _{k=1}^s}{\widehat{a}}_{ik}{H}_{Z_{G_j}^o}\left(y_{n_k},{\delta }_{n_k},{\tilde{p}}_{n_k}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}i=1,\cdots ,s\end{array}$(24)

The fifth equation of (20) can be written in the following explicit form:

$\begin{array}{l}{H}_{{\delta }_{x_j}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)=0,H_{{\delta }_{l_j}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)=0\\ H_{{\delta }_{m_j}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)=0,H_{{\delta }_{n_j}}\left(y_{n_i},{\delta }_{n_i},{\tilde{p}}_{n_i}\right)=0\end{array}$(25)

4.5. Fourth Order Conditions for Symplectic PRK Schemes

A fundamental tool in the study of order-$s$ conditions is the theory of B-series and the associated rooted-tree calculus, following Butcher [13]. For Partitioned Runge-Kutta schemes, this theory extends to P-series and the corresponding calculus on bicolored rooted trees. The latter allows the order-s conditions to be expressed in terms of the coefficients $a_{ik}$ , $b_i$ , ${\widehat{a}}_{ik}$ and ${\widehat{b}}_i$ . We assume the reader has basic background on B-series and P-series [15]. The derivation of the order-$s$ conditions for the symplectic Runge-Kutta method relies on graph theory, as shown by Julien Laurent-Varin and Hager [14] [28]. We will begin by deriving the order conditions for $s=4$ . For the table of order conditions for $s\ge 4$ , the reader may refer to the work of J. Laurent-Varin [14]. Defining $d_k={\displaystyle \sum }_{i=1}^s{b}_i{a}_{ik}$ and $c_i={\displaystyle \sum }_{k=1}^s{a}_{ik}$ , the 4th order conditions are collected in the following table:

Table 1. Equations of the eight 4th order conditions of SPRK algorithm [14].

$G_i$	Graph	Condition	$G_i$	Graph	Condition
$G_1$		${\displaystyle \sum }_{l,k=1}^4{a}_{lk}{d}_k{d}_l=\frac{1}{8}$	$G_2$
$G_3$		${\displaystyle \sum }_{i,k=1}^4\frac{b_i}{b_k}{a}_{ik}{c}_i{d}_k=\frac{5}{24}$	$G_4$		${\displaystyle \sum }_{k,i=1}^4{b}_i{a}_{ik}{c}_i{c}_k=\frac{1}{8}$
$G_5$		${\displaystyle \sum }_{k=1}^4{c}_k^2{d}_k=\frac{1}{12}$	$G_6$		${\displaystyle \sum }_{k=1}^4{b}_i^2{c}_i^3=\frac{1}{4}$
$G_7$		${\displaystyle \sum }_{k=1}^4\frac{1}{b_k}{c}_k{d}_k^2=\frac{1}{12}$	$G_8$		${\displaystyle \sum }_{l=1}^4\frac{1}{b_l^2}{d}_l^3=\frac{1}{4}$

4.6. Simulation Setup

The simulation is carried out over a fixed time horizon of 200 seconds, with a discretization step of 0.025. The initial state vector of the leader vehicle is set to $y_1\left(0\right)=\left(3000,22,2.5,22,2.5,1.5,0,0,0,5,5,5,50,20,0\right)$ , while initial state of the follower vehicle is set to $y_2\left(0\right)=\left(3250,22,2.59,22,3.5,1.5,0,0,0,5,5,5,125,21,0\right)$ . The initial control vector, shared by both vehicles, is defined as ${\delta }_0=\left(0.2,0,0.02,0\right)$ . Regarding constraint handling, the penalty coefficients for state-variable violations are sampled within the interval $\left(0,0.1\right)$ for all relevant state constraints, while the optimal penalty parameter $\rho$ is fixed at 1.5 × 10⁻³. For the objective function, the coupling parameters $\lambda$ are generated uniformly over $\left(0,1\right)$ . The penalty coefficients are fixed to ${\eta }_1=0.1$ , ${\eta }_2=0.025$ and ${\eta }_3=0.02$ respectively and the separation thresholds are set to $S_{\text{min}}=10$ m and $Y_{\text{max}}=1.5$ m. The numerical solver is configured with the following settings: a convergence tolerance (tol) of 10⁻⁴, a maximum number of iterations (max_iter) of 500, a constraint violation tolerance (constr_viol_tol) of 2.5 × 10⁻⁴, a dual infeasibility tolerance (dual_inf_tol) of 0.6 × 10⁻⁴, an adaptive barrier strategy (mu_strategy = adaptive), no NLP scaling (nlp_scaling_method = none), and error halting enabled (halt_on_ampl_error yes).

5. Digital Results of SPRK4 Method

The data and parameters of the model are extracted from the technical sheets of two different vehicles, the Land Cruiser 204ch-BVA6, 5 doors-5 seats, 2.8 L D-4D Diesel and the Toyota Hilux IV(2) X-TRA Cab 204ch-BVA6 4 doors-4 seats, 2.8 L D-4D Diesel, while other data are sourced from the literature [29]-[32]. All these data are presented in the Appendix (See Table A1 and Table A2).

Note that in this work, the inequality constraints are handled through a penalization approach directly embedded in the cost functional, which penalizes trajectories that violate the prescribed bounds on the state variables. Feasibility is assessed a posteriori by verifying that the computed trajectory remains within the admissible region at each discretization step, without explicitly enforcing complementarity conditions. It should be noted that, in our setting, only penalization is enforced, no complementarity condition of the form $\mu \cdot g\left(x\right)=0$ is explicitly imposed which simplifies the numerical resolution while still ensuring that constraint violations are sufficiently discouraged through the choice of the penalization parameter.

5.1. Numerical Results for Symplectic PRK4 Schemes

Taking into account the fourth order conditions from Table 1 for the practical implementation of the SPRK algorithm, the numerical computation of the coefficients used in the MATLAB code produces the following results:

Iteration	Func-count	Residual	First-Order optimality	Lambda	Norm of step
0	401	73511.3	3.28e+05	0.01	−
1	802	6687.11	1.04e+04	0.001	7.16992
2	1203	72.1339	746	0.0001	0.972471
3	1604	20.7787	8.13e+03	1e−05	0.357347
4	2007	2.94188	1.57e+03	0.001	1.71639
5	2409	1.39244	159	0.01	0.442513
6	2810	1.22495	1.05e+03	0.001	0.579999
7	3211	0.264347	383	0.0001	0.920228
8	3612	0.0186675	36.6	1e−05	0.333117
9	4013	0.000765176	4.33	1e−06	0.122321
10	4414	3.22118e−06	0.241	1e−07	0.0235201
11	4815	7.11017e−11	0.00113	1e−08	0.00158433
12	5216	1.69507e−20	1.66e−08	1e−09	7.43876e−06

$A=\left[\begin{array}{cccc}-1.139217& -3.376899& 7.154047& 0.468846\\ 1.608494& -3.662597& 1.611763& -0.030868\\ 0.523716& -3.335599& -0.903264& 1.910918\\ 1.494769& -6.200461& 2.325248& 1.551073\end{array}\right]$

$b=\left[\begin{array}{c}-0.000037\\ 0.004776\\ -0.000721\\ -0.005069\end{array}\right]$

$\widehat{A}=\left[\begin{array}{cccc}1.13918& 207.63092& -10.20611& -204.78842\\ -0.026198& 3.667373& -0.5042735& -6.585918\\ -0.367166& 10.68131& 0.902543& -16.352756\\ -0.003459& -0.024308& -0.272524& -1.556169\end{array}\right]$

5.2. Pontryagin’s Maximum Principle Outcomes

The following results were obtained using the indirect approach based on the Pontryagin’s Maximum Principle, implemented in the AMPL, “A Mathematical Programming Modeling Language”, environment and solved with the open-source solver IPOPT. The final feasibility and optimality error margin is: “EXIT= Optimal Solution Found, Dual infeasibility = 4.06e−006, Constraint violation = 2.52e−005, Complementarity = 1.005e−013, Overall NLP error = 2.52e−005”.

Figure 1 shows the evolution of the overall traffic flow metric $J_{G,12}$ , whose values range between 0.3 and 1.69, as well as the primary infeasibility errors Inf_pr and dual infeasibility errors Inf_du, which converge to 0. This result indicates a fluid traffic state in which vehicles travel normally at high speeds.

Fluidity function $J_{G,12}$ is measured on a scale from 0 to 1, occasionally exceeding 1 depending on the normalization adopted, where values below 0.5 indicate dense, congestion-prone traffic likely to generate severe delays, while values above 0.5 reflect a transition toward a smoother flow regime. In this case, the evolution of $J_{G,12}$ reveals that traffic is progressing toward a free-flow state at high speeds, indicating that queues are well controlled and dissipate without significant buildup. In practice, this result suggests that intervention efforts should be strategically concentrated on low-fluidity zones, where targeted actions can be deployed rapidly and effectively to prevent congestion from developing.

The results in Figure 2 show that the assumption of variation in the vehicle masses $m_1$ and $m_2$ remains valid, which supports an integrated approach adapted to contemporary transportation challenges. The aerodynamic speed $v_{a1}$ varies from 22 m/s to 3.75 m/s² and $v_{a2}$ decreases from 22 m/s to 21.89 m/s. The aerodynamic acceleration ${\gamma }_{a1}$ varies from 2.5 m/s² to 2.486 m/s² and ${\gamma }_{a2}$ varies from 2.59 m/s to 2.576 m/s². These results emphasize the importance of avoiding collision risks due to interactions between vehicles on the road. This may be desirable for several reasons: maintaining passenger comfort and avoiding sudden maneuvers; promoting cautious or fuel-efficient driving aimed at minimizing fuel consumption and reducing vehicle wear; serving as a precautionary indicator to prevent damage; and reducing the risk of losing control of the vehicle.

Figure 1. Evolution of the objective function and infeasibility errors.

Figure 2. Evolution of masses, the aerodynamic velocities, and accelerations of the vehicles.

Figure 3 shows that the primary thrust control ${\delta }_x$ varies from 0.1 to 0.8 for both vehicles and also exhibits a steady behavior around the value 0.8. The roll control ${\delta }_{l1}$ ranges from −8.48142 × 10⁻¹⁵ to 0.00321695, while ${\delta }_{l2}$ ranges from −2.4446 × 10⁻¹⁶ à −4.11088 × 10⁻⁸, which demonstrates the stability of the second vehicle and indicates safe, secure driving maneuvers. The pitch control ${\delta }_{m1}$ varies from 0 to 0.02, whereas ${\delta }_{m2}$ varies from 0 to 0.095. The yaw command ${\delta }_{n1}$ ranges from 0 to 0.35, while ${\delta }_{n2}$ ranges from 0 to −0.036. The results obtained for pitch and yaw controls demonstrate directional stability for both vehicles and suggest stable and safe driving, with a good balance of forces and effective vehicle handling, which is essential for safety and comfort.

Figure 3. Controls of the vehicles.

Figure 4 shows that the longitudinal speed $u_1$ decreases from 22 m/s to 21.8603 m/s while $u_2$ varies from 22 m/s to 25.7809 m/s. The lateral speed $v_1$ varies from 3 m/s to 5 m/s whereas $v_2$ decreases from 3.5 m/s to 3.49678 m/s. The vertical speeds $w_1$ and $w_2$ both slightly decrease around 1.5 m/s to approximately 1.49221 m/s for the two vehicles. The roll rate $p_1$ ranges from 0/s to 0.0332393/s while $p_2$ ranges from 0/s to 0.000424323/s. According to SAE and ISO standards, these results indicate upstream congestion slots and a situation where the vehicles experience an increase in roll associated with a lateral maneuver.

Figure 4. Time evolution of the aerodynamic velocity components and the roll rate of the vehicles.

Figure 5 shows that the pitch rate $q_1$ evolves from −5.45 × 10⁻⁷˚/s to 3.37 × 10⁻⁷ while $q_2$ varies from 0˚/s to −5.48 × 10⁻⁹˚/s. The yaw rate $r_1$ ranges from 0˚/s to 3.427 × 10⁻⁷˚/s while $r_2$ ranges from 0˚/s to 1.01 × 10⁻⁷˚/s. The roll angles ${\theta }_1$ and ${\theta }_2$ undergo a slight decrease from 5˚ à 4.97˚ for both vehicles. The pitch angles ${\varphi }_1$ vary from 5˚ to 5.11˚ while ${\varphi }_2$ vary from 5˚ to 4.97˚. These results reveal a decrease in pitch accompanied by an increase in yaw, indicating an interaction between longitudinal dynamics and lateral maneuvers. This suggests that the vehicles reduce their lateral inclination while adjusting their nose-down/up attitude to follow a new trajectory, requiring an analysis of control actions for precise diagnosis.

Figure 5. Time evolution of pitch rate, yaw rate, roll angle, and pitch angle.

Figure 6 shows that the yaw angle $\psi$ varies slightly from 5 to 4.97 for both vehicles. The longitudinal and lateral separations are 75 m and 1 m, respectively. The position $X_{G,1}^o$ increases from 50 m to 58.14 m while $X_{G,2}^o$ varies from 125 m to 133.5 m for the second vehicle. The position $Y_{G,1}^o$ changes from 20 m to 20.1149 m while $Y_{G,2}^o$ varies from 21 m to 21.1069 m for the second vehicle. The position $Z_G^o$ varies from 0m to 0.0058 m for both vehicles. These results show a very small negative change in the yaw angle, reflecting a very slight negative-directed rotation of the vehicle. The position variations indicate that the vehicles move forward relative to an observer attached to the inertial frame.

Figure 6. Yaw angle evolution and 3D position of the centers of mass G1 and G2.

Figure 7 shows the evolution of the components of the adjoint state vector. The adjoint state is an additional dimensionless variable that can be interpreted as a dynamic multiplier measuring the marginal value of the state with respect to the cost. Some corresponding components are constant for both vehicles and oscillate only around the values −1, 0, and 1, this result indicates numerical stability.

Figure 7. Evolution of the adjoint state components.

To further highlight the effectiveness of the penalized PMP approach, a comparison with a standard non-symplectic discretization applied to the same scenario reveals a clear advantage in terms of solution quality and dynamical consistency. Unlike conventional discretization schemes, which do not inherently respect the geometric structure of the optimal control problem and may thus yield suboptimal or numerically inconsistent trajectories, the penalized PMP approach is rigorously grounded in Pontryagin’s necessary conditions of optimality. Although this theoretical rigor comes at the cost of increased implementation complexity requiring the resolution of a two-point boundary value problem and careful handling of the adjoint dynamics. It ultimately delivers a more reliable and physically meaningful control law, particularly in the presence of state constraints enforced through penalization.

6. Conclusions

This paper demonstrates that Pontryagin’s Maximum Principle provides a powerful foundation for optimizing the dynamics of vehicle queues in road traffic systems. By integrating traffic-flow modeling with optimal control theory, we derive optimal feedback strategies under a penalized joint formulation that significantly reduce delays and enhance throughput while accounting for constraints and environmental factors.

Beyond the conceptual framework, two ingredients prove especially important for practical and reliable performance. First, exterior penalization plays a crucial role in enforcing the feasibility and robustness of the optimized queue dynamics, ensuring that the numerical solution remains physically meaningful and constraint-compliant. Second, the use of symplectic PRK schemes is key to obtaining accurate and stable long-horizon simulations of the resulting Hamiltonian (state-adjoint) system. In particular, these structure-preserving integrators improve numerical behavior and reliability when computing optimal strategies for complex traffic scenarios.

Overall, the combination of PMP-based indirect optimization with exterior penalization and symplectic PRK discretizations offers an effective methodology for advancing queue management and traffic efficiency, and it provides useful insights for future developments in transportation optimization, extension to larger networks, tighter constraint handling, and real-time implementations.

Acknowledgements

The first author acknowledges that this research received no external funding.

Appendix

Table A1. Model parameter values.

		Value
Denomination	Notation	Vehicule 1	Vehicule 2
Vehicle inertia moments (kgm−2)	$I_{xx}$	1739.3125	1814.2
	$I_{yy}$	6707.41	8581.35
	$I_{zz}$	6744.71	8621
Vehicule length	$L$	4.84 m	5.33 m
Vehicule width	$t$	1.885 m	1.85 m
Vehicule height	$H$	1.845 m	1.81 m
Motor mechanical efficiency	${\eta }_m$	≈0.9	≈0.9
Fuel Conversion efficiency	${\eta }_c$	≈0.4	≈0.4
Basic aerodynamic drag coefficient	$C_{a_0}$	≈0.35	≈0.35
Basic lateral force coefficient	$C_{y_0}$	≈0.1	≈0.1
Basic vertical force coefficient	$C_{z_0}$	≈-0.1	≈-0.1
Basic Roll moment coefficient	$C_{l_0}$	≈0.1	≈0.1
Fuel cosumption coefficient	$C_{cse}$	≈1 × 10−4 kg/m	≈1 × 10−4 kg/m
Internal motor power	$P_i$	2 × 105 W	2.05 × 105 W
Throttle ejection speed	$v_e$	30 m/s	30.5 m/s
Calorific value of fuel	$LH{V}_{fuel}$	≈42.7 × 105 J/kg	≈42.7 × 105 J/kg
Gas outlet pressure	$P_e$	1.0 × 105 Pa	1.01 × 105 Pa
Exhaust outlet surface	$A_e$	≈0.1 m2	≈0.12 m2
Vehicle frontal area	$A_f$	2.89 m2	3.01 m2

Table A2. Limit values for Vehicle state variables and controls.

Denomination	Minimum value	Maximum value
Mass of the vehicle 1	$m_{10}=3000\text{kg}$	$m_{1f}=2500\text{kg}$
Mass of the vehicle 2	$m_{20}=3250\text{kg}$	$m_{2f}=2500\text{kg}$
Aerodynamc vehicle 1 speed	$v_{a_1\min }=0\text{m}/\text{s}$	$v_{a_1\max }=55\text{m}/\text{s}$
Aerodynamc vehicle 2 speed	$v_{a_2\min }=0\text{m}/\text{s}$	$v_{a_2\max }=50\text{m}/\text{s}$
Aerodynamc vehicles acceleration	${\gamma }_{a_1\min }={\gamma }_{a_2\max }=-5\text{m}/{\text{s}}^2$	${\gamma }_{a_1\max }={\gamma }_{a_2\max }=5\text{m}/{\text{s}}^2$
Longitudinal vehicle 1 speed	$u_{1\min }=22\text{m}/\text{s}$	$u_{1\max }=55\text{m}/\text{s}$
Longitudinal vehicle 2 speed	$u_{2\min }=22\text{m}/\text{s}$	$u_{2\max }=45\text{m}/\text{s}$
Lateral vehicles speed	$v_{1\min }=v_{2\min }=-5\text{m}/\text{s}$	$v_{1\max }=v_{2\max }=5\text{m}/\text{s}$
Vertical vehicles speed	$w_{1\min }=w_{2\max }=-3\text{m}/\text{s}$	$w_{1\min }=w_{2\max }=3\text{m}/\text{s}$
Vehicle Roll velocity	$p_{1\min }=p_{2\min }=-5˚/\text{s}$	$p_{1\max }=p_{2\max }=5˚/\text{s}$
Vehicle Pitch velocity	$q_{1\min }=q_{2\min }=-2˚/\text{s}$	$q_{1\max }=q_{2\max }=2˚/\text{s}$
Vehicle Yaw velocityt	$r_{1\min }=r_{2\min }=-2˚/\text{s}$	$r_{1\max }=r_{2\max }=2˚/\text{s}$
Roll angles	${\theta }_{1\min }={\theta }_{2\min }=-15˚$	${\theta }_{1\max }={\theta }_{2\max }=15˚$
Pitch angles	${\varphi }_{1\min }={\varphi }_{2\min }=-3˚$	${\varphi }_{1\max }={\varphi }_{2\max }=30˚$
Yaw angles	${\psi }_{1\min }={\psi }_{2\min }=0˚$	${\psi }_{1\max }={\psi }_{2\max }=45˚$
Vehicle longitudinal position	$X_{G_1\min }=X_{G_2\min }=0\text{m}$	$X_{G_1\max }=X_{G_2\max }={10}^3\text{m}$
Vehicle lateral position	$Y_{G_1\min }=Y_{G_2\min }=-{10}^2\text{m}$	$Y_{G_1\max }=Y_{G_2\max }={10}^2\text{m}$
Vehicle vertical position	$Z_{G_1\min }=Z_{G_2\min }=0\text{m}$	$Z_{G_1\max }=Z_{G_2\max }=10\text{m}$
Thrust control	${\delta }_{x_1\min }={\delta }_{x_2\min }=0.1$	${\delta }_{x_1\max }={\delta }_{x_2\max }=0.8$
Roll control	${\delta }_{l_1\min }={\delta }_{l_2\min }=-0.0275$	${\delta }_{l_1\max }={\delta }_{l_2\max }=0.0275$
Pitch control	${\delta }_{m_1\min }={\delta }_{m_2\min }=0$	${\delta }_{m_1\max }={\delta }_{m_2\max }=0.095$
Yaw control	${\delta }_{n_1\min }={\delta }_{n_2\min }=-0.036$	${\delta }_{n_1\max }={\delta }_{n_2\max }=0.35$