Optimization Algorithms for Distributed Model Predictive Control: Benchmarking Dual Decomposition Approaches

This paper introduces a benchmark comparison of distributed optimization algorithms based on dual decomposition for model predictive control problems with constraint coupling. Such problems can be viewed as multiple subsystems linked through constraints on shared, limited resources. In a dual decomposition approach, the allocation and use of these resources are coordinated by repeatedly calculating their associated prices and distributing those values to the participating subsystems. After a short overview of model predictive control, various architectures and communication structures suitable for distributed implementations are outlined. The centralized control problem is then decomposed into several subproblems by introducing dual variables. A range of distributed optimization methods is examined, including the recently developed quasi-Newton dual ascent algorithm. Additionally, an epigraph reformulation of the bundle, combined with a line search technique, is proposed for the quasi-Newton dual ascent method; these enhancements improve its numerical stability and accelerate convergence relative to the earlier trust-region approach. Lastly, the quasi-Newton dual ascent algorithm is benchmarked against the subgradient method, the bundle trust method, and the alternating direction method of multipliers across a broad set of test problems. All benchmark instances are openly accessible on GitHub. 

The growing complexity of industrial systems has driven the adoption of sophisticated control techniques. For many years, decentralized control structures relying on numerous independent PID loops have been widely implemented in practice [1]. These approaches typically presume minimal or no interaction among the control loops within the overall system. Nevertheless, the pursuit of greater efficiency in large-scale systems has resulted in stronger interconnections among components, thereby weakening the effectiveness of purely decentralized strategies. Moreover, contemporary systems—and their corresponding control systems—are required to meet a variety of operational and safety constraints while simultaneously optimizing multiple performance objectives. Model predictive control (MPC) is an optimization-based strategy that uses a plant model to forecast future system behavior from the current state and prospective control actions. Formulating the control task as an optimization problem enables direct minimization or maximization of performance criteria while explicitly enforcing all relevant constraints. Today, problems involving hundreds or even thousands of variables can be solved efficiently, largely thanks to major advances in the performance of optimization solvers. Despite these developments, certain scenarios still call for more decentralized control strategies. A prominent example involves large-scale systems composed of multiple autonomous subsystems, each having its own goals and limitations. Although these subsystems can make independent decisions, they may remain coupled—for instance, through shared energy or material flows. In such cases, global constraints must be upheld even when local control problems are addressed in a distributed fashion. One key motivation for this architecture is the need to preserve confidentiality among subsystems, particularly when they belong to separate organizations unwilling or unable to share details of their internal dynamics and local constraints. Distributed optimization solves the problem of coordinating individual subsystem controllers, leading to distributed model predictive control (DMPC) problems.

A class of distributed optimization algorithms particularly well-suited for coordinating resource-coupled DMPC problems is dual-decomposition-based methods. These approaches are frequently called price-based coordination schemes, in which the introduced dual variables are interpreted as prices associated with shared, limited resources; these prices guide the subsystems toward solutions that satisfy the overall system-wide constraints. 

The present paper makes two main contributions: First, the recently introduced quasi-Newton dual ascent algorithm [2, 3] is enhanced by reformulating its cutting planes as linear constraints through an epigraph representation. In addition, a line search procedure with constraint handling is incorporated, which accelerates convergence, especially near the optimal solution. Second, a comprehensive set of DMPC benchmark problem instances has been created and released for public use on GitHub. Although several prior studies have benchmarked DMPC algorithms, many are limited to single instances [4-6] or do not provide open access to the benchmark problems [7]. 

In this work, vectors are denoted by bold lowercase letters (𝐱) and matrices by bold uppercase letters (𝐀). In model predictive control, time is typically discretized into uniform sampling intervals of length 𝑇𝑠, resulting in discrete time points 𝑡𝑘 = 𝑡0 + 𝑘 ⋅ 𝑇𝑠,, where 𝑡0 denotes the current time. The symbols 𝐱𝑘 and 𝐮𝑘 denote the plant states and control inputs, respectively, at time 𝑡𝑘. The value of a variable 𝐳 at iteration 𝑡 of a dual decomposition-based algorithm is indicated by (𝑡). The Euclidean norm is written as ‖ ⋅ ‖₂. 

MPC methods generally rely on a discrete-time model to forecast the future evolution of the plant states 𝐱 ∈ ℝ^{𝑛ₓ} based on the current states and control inputs 𝐮 ∈ ℝ^{𝑛ᵤ}. A standard MPC optimization problem is formulated as shown in Eq. 1. 

Constraint (Eq. 1b) describes a discrete-time dynamic model that predicts the plant state 𝐱𝑘₊₁ at the subsequent sampling instant from the current state 𝐱𝑘 and applied input 𝐮𝑘, with system matrices 𝐀 ∈ ℝ^{𝑛ₓ×𝑛ₓ} and 𝐁 ∈ ℝ^{𝑛ₓ×𝑛ᵤ}. Predictions of future states and inputs are generated over a prediction horizon of length 𝑁ₚ. It should be noted that control inputs need not be optimized over the full prediction horizon, particularly when computational cost is a concern. In such cases, inputs may be optimized only over a shorter control horizon 𝑁𝑐 < 𝑁ₚ, with the final computed input held constant for the remaining steps (i.e., 𝐮𝑁𝑐 = 𝐮𝑁𝑐₊₁ = ⋯ = 𝐮𝑁ₚ₋₁). For this study, it is assumed that 𝑁ₚ = 𝑁𝑐. Constraint (Eq. 1c) specifies the initial conditions of the plant, where 𝐱̃(𝑡₀) represents the measured states at the current time 𝑡₀. In situations where not all states are directly measurable, state estimation techniques may be employed to reconstruct them from available outputs; however, state estimation lies beyond the scope of this paper. An overview of state estimation techniques in distributed frameworks is provided by Necoara et al. [8]. 

One of the primary strengths of MPC over alternative control strategies is its ability to explicitly incorporate constraints on both states and inputs directly into the optimization formulation [9]. These constraints are collectively expressed in (Eq. 1d) for the states and in (Eq. 1e) for the inputs.

At every sampling instant, the controller calculates a sequence of future inputs across the prediction horizon 𝑁𝑝 by solving the open-loop optimization problem defined in (Eq. 1). 

A frequent goal in model predictive control is to make the system follow a specified reference signal throughout the horizon while keeping control actions as moderate as possible, which leads to the cost function.

Where 𝐇𝐱 ∈ ℝ^{𝑛ₓ×𝑛ₓ} and 𝐇𝐮 ∈ ℝ^{𝑛ᵤ×𝑛ᵤ} are symmetric positive (semi-)definite weighting matrices. The expression in Eq. 2 is usually called the stage cost, in contrast to the terminal cost (𝐱𝑁ₚ). One standard form of the terminal penalty penalizes the squared distance between the states at the final time step and their desired reference values, namely. 

Careful tuning of the terminal cost term is essential to ensure satisfactory closed-loop behavior and stability [9]. Discussion of these design issues, however, is beyond the objectives of the current study. A key strength of MPC lies in its ability to handle more advanced control objectives than mere setpoint tracking by embedding them directly into the optimization setup [10]. For example, one can optimize process economics directly, leading to the field known as economic model predictive control (EMPC). For a thorough treatment of EMPC, the reader is referred to the review by Ellis et al. [11]. 

After the open-loop optimization has been completed, only the initial control move 𝐮₀ is actually implemented. Then, at the next sampling time, the plant’s new state is measured (or estimated), and the optimization cycle restarts.

Because MPC requires an optimization problem to be solved repeatedly at each discrete time step, difficulties arise when that problem becomes excessively large—either exceeding available computation time or preventing critical data from being gathered at a single central location due to a multi-subsystem architecture. One effective way to overcome these limitations is to break the global MPC task into smaller pieces and solve them in a distributed way.

Research on distributed MPC (DMPC) has primarily focused on two distinct problem classes.

The first class involves a collection of 𝑁ₛ dynamically interconnected subsystems indexed by the set 𝐼 = {1, …, 𝑁ₛ}. The corresponding optimization problems take the form 

With 𝐱ᵢ and 𝐮ᵢ standing for the state and input vectors of the 𝑖-th subsystem. The complete vectors of all subsystem states and inputs are assembled as 𝐱 = [𝐱₁ᵀ … 𝐱_{𝑁ₛ}ᵀ]ᵀ and 𝐮 = [𝐮₁ᵀ … 𝐮_{𝑁ₛ}ᵀ]ᵀ. Each subsystem is equipped with its private terminal cost 𝐽ᵢ(𝐱ᵢ_{𝑁ₚ}) and per-step running cost 𝐽(𝐱ᵢₖ, 𝐮ᵢₖ). The matrices 𝐀_{𝑖𝑗} ∈ ℝ^{𝑛_{𝑥ᵢ}×𝑛_{𝑥ⱼ}} and 𝐁_{𝑖𝑗} ∈ ℝ^{𝑛_{𝑥ᵢ}×𝑛_{𝑢ⱼ}} quantify how the states and inputs of subsystem 𝑗 affect the evolution of subsystem 𝑖.

The second class features subsystems that interact solely via inequality constraints. 

In this setting, the individual subsystem dynamics (Eq. 5b) are completely decoupled. The only source of coupling is the collection of global resource constraints (Eq. 5f), which enforce limits on shared scarce resources. Each matrix 𝐑ᵢ ∈ ℝ^{𝑛ᵣ×𝑛_{𝑢ᵢ}} converts the inputs of subsystem 𝑖 into the net consumption or generation of the various resources, while the vector 𝐫 ∈ ℝ^{𝑛ᵣ} defines the upper bound on resource availability. Within these constraints, a positive entry indicates resource use, and a negative entry indicates resource supply.

The focus of the present article is exclusively on the second category—DMPC problems of type (5)—which lend themselves naturally to resolution through dual decomposition techniques in a distributed optimization setting. For completeness, it is worth noting that problems in the first category (4) can also be tackled using dual decomposition methods.

DMPC architectures 

Besides the specific form of interaction among the subsystems, the selected distributed model predictive control (DMPC) strategy is largely determined by the permitted communication constraints. Several DMPC schemes have been examined in detail by authors such as Scattolini [12] and Christofides et al. [13]; several representative examples appear in Figure 1. 

To keep the presentation straightforward, we temporarily neglect any influence that one subsystem’s control output exerts on others, which corresponds to setting 𝐁𝑖𝑗 = 𝟎 whenever 𝑖 ≠ 𝑗 ( Eq. 4b). 

The layout in Figure 1(a) represents the case in which DMPC controllers completely disregard interconnections among subsystems and therefore do not exchange any data. This arrangement is commonly referred to as decentralized MPC. Here, any coupling effects are treated as unknown disturbances that each local controller must reject on its own. Strong inter-subsystem coupling frequently leads to degraded closed-loop behavior under this scheme. In addition, purely decentralized DMPC cannot reliably enforce globally defined constraints because there is no mechanism for coordinating constraint satisfaction across subsystems. Further discussion of decentralized MPC approaches can be found in Christofides et al. [13]. 

Whenever the subsystems are coupled in a way that requires explicit handling inside the local optimizers, at least limited information sharing is essential. Figure 1(b) illustrates the scenario that allows unrestricted, direct communication links among all DMPC controllers. Through these links, the controllers can transmit predicted state and input trajectories, cost function values, active constraints, and other relevant quantities. 

Algorithms that adopt this fully connected communication pattern are commonly distinguished by the number of information-exchange steps executed during each MPC time step—i.e., the count of negotiation cycles before the computed inputs are actually implemented. Both single-round and multi-round communication strategies are systematically reviewed by Camponogara et al. [14]. 

The main practical consideration when choosing between one or several communication rounds per sampling interval is the CPU time needed to solve a single local DMPC problem. Multi-round methods—often described as negotiation among agents [15]—generally yield improved overall performance and, in many cases, converge toward the globally optimal solution that a centralized MPC would produce. The drawback is a significantly higher computational load, since each subsystem must resolve its optimization multiple times within a single control cycle. For applications that demand very short sampling periods, repeated negotiations become infeasible. A typical example is collision-free control of robotic arms, where the forecasted state sequence is usually forwarded to adjacent controllers [16]. 

In contrast, the configuration shown in Figure 1(c) forbids any peer-to-peer communication between the local DMPC controllers. Instead, all local optimizers report to and receive directives from a single coordinating entity, resulting in a hierarchical architecture. The coordinator gathers the local solutions, computes appropriate price (or Lagrange multiplier) updates, and broadcasts these signals back to the subsystems. Dual decomposition techniques rely precisely on this kind of information flow, with the coordinator transmitting updated dual variables at every iteration. 

Hierarchical DMPC proves especially advantageous for systems whose constraints couple multiple subsystems yet where direct controller-to-controller data sharing is restricted or prohibited—for reasons such as data privacy or organizational boundaries. One important characteristic of hierarchical schemes is that they invariably involve several rounds of communication between the local units and the coordinator during each sampling interval. For this reason, the method is realistic only when the plant exhibits slow dynamics and the sampling time is long enough to accommodate repeated solution of the optimization subproblems until convergence or feasibility is reached. Representative application domains include large-scale chemical processes and electric power networks. 

The subsequent section briefly surveys selected recent contributions to DMPC, focusing primarily on cooperative and hierarchical formulations.

 DMPC literature review 

Research into distributed control strategies for large-scale processes has remained vibrant since the 1970s [17]. In-depth surveys of distributed model predictive control (DMPC) appear in works by, among others, Scattolini [12], Necoara et al. [8], Christofides et al. [13], Negenborn and Maestre [18], and Müller and Allgöwer [19]. 

Stewart et al. introduce a cooperative DMPC technique that functions as an approximate centralized MPC and thereby establishes closed-loop stability [20]. Maestre et al. describe a cooperative DMPC method tailored to two agents involving two information-exchange steps per control interval [21]. The same group also develops a cooperative DMPC framework for multiple agents, in which, during each communication cycle, a randomly chosen subset of agents proposes control moves that neighboring agents adopt only if those moves improve the global performance index [15]. Zheng et al. outline a DMPC procedure designed for networks of interconnected systems featuring direct controller-to-controller links [22]. 

Approaches based on dual decomposition have addressed systems with coupling in both the dynamics and the constraints. Giselsson and Rantzer apply dual decomposition to dynamically coupled systems and introduce a termination condition that ensures limited suboptimality together with asymptotic stability [23]. In related work, Giselsson et al. propose an accelerated gradient technique for dual decomposition DMPC that incorporates L₁-norm regularization in the cost [24]. Giselsson and Rantzer further describe an early-stopping rule based on adaptive constraint relaxation, enabling faster convergence of dual decomposition iterations in dynamically coupled setups [25]. Köhler et al. investigate conditions for recursive feasibility when halting a dual decomposition distributed optimizer ahead of full convergence [26]. Doan et al. utilize dual decomposition combined with primal averaging and constraint adjustment to secure convergence in DMPC scenarios involving coupled dynamics and constraints alike [27]. 

Biegel et al. demonstrate the application of MPC to congestion control within an electricity network populated by smart consumers equipped with storage units [28]. This concept is later adapted to a fully distributed DMPC framework [29]. Biegel et al. also illustrate the broad applicability of dual decomposition DMPC to systems linked solely through constraints [30]. Pflaum et al. compare primal versus dual decomposition DMPC strategies for distributed energy management in smart neighborhoods [31]. In a separate study [32], they assess how well dual decomposition DMPC scales—using a bundle-based solver—when varying the count of subsystems, the quantity of constraints, and the largest bundle dimension. 

Razzanelli et al. implement DMPC with a subgradient algorithm to orchestrate energy exchanges across a cluster of microgrids [33]. Eser et al. adopt the alternating direction method of multipliers (ADMM) for decentralized regulation of building climate and energy systems [34]. Maxeiner and Engell employ ADMM to handle DMPC of semi-batch chemical reactors that share scarce input resources [35]. That effort is broadened in the work by Maxeiner and Engell [5], where the authors conduct a thorough benchmark of the subgradient approach, ADMM, and the augmented Lagrangian alternating direction inexact Newton (ALADIN) method for resource-constrained coupled problems. Houska and Shi offer an instructional overview of applying ALADIN to DMPC tasks [36]. Yfantis et al. introduce a quadratic model-based solver for DMPC under shared resource limits and compare its efficiency with the subgradient technique [37]. 

Conte et al. explore numerical properties of dual decomposition DMPC when accelerated gradient or ADMM solvers are used [7]. Their analysis specifically quantifies how factors such as coupling intensity, system stability, starting conditions, network scale, and topology affect runtime performance. Stomberg et al. benchmark six different DMPC solvers for dynamically interacting subsystems: the basic subgradient method, an accelerated subgradient variant, ADMM, a distributed active-set procedure, a nearly decentralized interior-point algorithm, and Jacobi-style iterations [6]. 

Beyond ADMM, DMPC also uses additional proximal splitting methods. Halvgaard et al. apply the Douglas-Rachford splitting scheme to coordinate distributed smart energy infrastructures [38]. A comprehensive compilation of DMPC methods, along with practical case studies, is available in the textbook edited by Maestre and Negenborn [39]. 

The present section focuses on DMPC grounded in dual decomposition principles. By decomposing an overarching constraint-coupled MPC formulation through the introduction of dual (Lagrange) multipliers, one obtains a dual optimization problem amenable to a fully distributed solution. Various iterative solvers for this dual problem are outlined, among them the quasi-Newton dual ascent procedure. For completeness, the section also illustrates how dual decomposition extends naturally to DMPC problems featuring dynamic coupling between subsystems.

Dual decomposition of constraint-coupled DMPC problems 

This section explains how dual decomposition is applied to separate the overall constraint-coupled model predictive control (MPC) problem into independent subproblems. By associating dual variables 𝝀𝑘 ∈ R𝑛𝐫 with each sampling instant 𝑘, the Lagrangian of the centralized problem can be expressed as follows. 

The full Lagrangian then splits additively into a sum of local Lagrangian terms, one for each subsystem. 

As a result, the dual function 

Becomes computable in a fully distributed fashion by independently solving the local DMPC optimization problems 

for every subsystem 𝑖. 

The core idea of dual decomposition-based optimization is to tackle the dual problem 

This naturally lends itself to parallel and distributed implementation, thanks to the dual function’s separability. The equality constraint (Eq. 10b) arises directly from the Karush–Kuhn–Tucker (KKT) stationarity conditions applied to the inequality-coupled constraints [40]. 

Several key theoretical features characterize both the dual function and the associated dual optimization problem. First, the optimal value attained by the dual problem (10) invariably serves as a lower bound to the optimal value of the original (primal) problem (5)—a relation known as weak duality. When the primal problem is convex and satisfies appropriate constraint qualifications, the primal and dual optimal values coincide exactly; this is termed strong duality [40]. Moreover, once the optimal dual variables are known, solving the local subproblems (9) at those dual values yields a globally optimal solution to the primal problem. 

Second, irrespective of the convexity properties of the primal problem, the dual function is always concave. That said, it is generally nonsmooth: the gradient with respect to the dual variables may fail to exist at certain points [3]. Consequently, the dual problem (10) constitutes a nonsmooth convex optimization task. The next section surveys numerical methods specifically designed to solve such nonsmooth dual problems. 

 Dual decomposition-based algorithms 

The current section describes the distributed optimization techniques based on dual decomposition that form the backbone of the methods presented in this work. It begins with a short overview of the classical subgradient method, followed by the bundle trust-region approach and the alternating direction method of multipliers (ADMM). Particular attention is then devoted to the more recent quasi-Newton dual ascent algorithm introduced by Yfantis et al. [2, 3]. Compared with the version detailed by Yfantis et al. [3], the algorithm described here incorporates two important enhancements: a reformulation of the inner update subproblem using an epigraph representation and the addition of a line-search procedure applied in the neighborhood of the solution.

 Subgradient method 

The most straightforward technique for addressing the nonsmooth dual optimization problem is the subgradient method [41]. This approach extends the notion of a gradient to functions that are not differentiable everywhere. For convenience, we represent the complete vector of dual variables across the entire prediction horizon as a single column vector. 

A subgradient of the dual function (Eq. 8) at iteration 𝑡 is obtained by evaluating the coupling constraints on the shared limited resources (Eq. 5f), yielding: 

Where 𝐮𝑘, (𝑡+1) refers to the control sequences determined by each subsystem through solution of the local problem (9) during iteration 𝑡 of the distributed procedure. 

In the subgradient method, the dual variables are then adjusted at every step by moving in the direction of this subgradient, according to 

Where [⋅]+ indicates projection onto the nonnegative orthant. The step-size coefficient (𝑡) is typically adjusted dynamically as iterations progress. The procedure stops once both the primal residual norm 

And the dual residual norm 

Falls below a user-specified tolerance, indicating that the original primal problem has become feasible and the dual variables have stabilized at a stationary point. 

Bundle trust method 

Unlike the subgradient method, which relies solely on the most recent subgradient to update the dual variables, bundle methods form a broader, often more effective family of solvers for nonsmooth optimization problems [42]. These algorithms accumulate subgradient information from several past iterations to guide the updates. 

Specifically, in each iteration, the following information is retained: 

Where 𝑛𝝀 = 𝑛𝐫 ⋅ (𝑁𝑝 − 1) gives the total dimension of the dual variable vector. The set (𝑡) is called the bundle; it comprises the dual variable values, corresponding subgradients, and dual function evaluations collected from earlier steps. To prevent excessive memory usage, only the most recent 𝜏 iterations’ data are kept in the bundle. 

Using this accumulated bundle, a piecewise linear outer approximation of the dual function—known as a cutting-plane model—is constructed and then maximized. Because this model remains nonsmooth, its optimization is reformulated as an epigraph problem, and a search direction for the dual variables is computed under an additional trust-region constraint. This results in the bundle trust method (BTM): 

With the linearization error defined as 

Once a promising direction has been determined, the dual variables are updated via 

Alternating direction method of multipliers 

ADMM introduces quadratic penalty terms that effectively smooth the nonsmooth dual function and promote convexity in the response surface [43]. When multiple subsystems are interconnected solely through shared, scarce resources, the optimal-exchange variant of ADMM is particularly suitable [44]. 

Under this formulation, auxiliary variables 𝐳𝑖𝑘 are introduced for every subsystem and time instant, allowing the construction of an augmented Lagrangian as follows: 

Theoretical analysis confirms that ADMM converges reliably when the penalty parameter 𝜌 is held constant. Nevertheless, practical implementations frequently achieve faster and more stable progress by adjusting 𝜌 dynamically across iterations. The adaptation scheme implemented here follows the strategy outlined by He et al. [45] and Wang and Liao [46]: 

The constants, 𝜏^incr and 𝜏^decr (each greater than 1), act as tunable hyperparameters. Once this rule is applied, the algorithm proceeds by sequentially updating the primal variables, dual variables, and auxiliary variables: 

One characteristic advantage of ADMM is that the augmentation term actively encourages the subsystem solutions to satisfy the primal constraints. In problems with inequality coupling, however, the algorithm may reach a point at which the augmented constraints are satisfied while a residual penalty contribution persists. At true optimality, this residual should disappear, so that the final dual variables produce a feasible point even when evaluated with the original (unaugmented) Lagrangian. 

To overcome this limitation, Maxeiner [47] recommended replacing the standard auxiliary update step (Eq. 22b) with a dedicated optimization subproblem: 

The consequences of this change are visualized in Figure 2. Employing the original update rule (Eq. 22b) typically causes ADMM to settle at dual variables that enforce resource limits only when the penalty is active; removing the augmentation then reveals violations of the shared constraints. This behavior appears in the left panel of Figure 2 for a pair of representative resources. 

In contrast, applying the proposed update (Eq. 23) drives ADMM toward the genuine dual optimum, delivering dual variables that guarantee feasibility in both the augmented and the unaugmented settings. The right panel of Figure 2 demonstrates this superior outcome. Consequently, all further discussion and results in this work rely on updating the auxiliary variables via (Eq. 23). 

Quasi-Newton dual ascent 

Yfantis et al. recently introduced the quasi-Newton dual ascent (QNDA) method [2, 3]. The approach relies on constructing a quadratic model to approximate the concave dual function. 

This strategy draws inspiration from classical Newton-type methods. Since the dual function is nonsmooth, however, exact gradients and Hessians cannot be computed. The algorithm therefore employs subgradients and updates an approximate Hessian at each step using the Broyden–Goldfarb–Shanno–Fletcher (BFGS) formula: 

Where, 

Denotes the change in the dual variables between consecutive iterations and: 

Represents the corresponding change in the subgradients. Once the quadratic model has been updated, it is maximized to obtain the next iterate of the dual variables. 

Although the model  provides a smooth quadratic surrogate, the true dual function remains nonsmooth, leading to substantial modeling inaccuracies—particularly near or at points where differentiability fails. To mitigate this drawback, the algorithm incorporates information from the bundle defined in (Eq. 16). While bundle trust methods (BTM) exploit the bundle to form a piecewise linear outer approximation of the dual function, the QNDA scheme instead uses the collected data to generate supporting hyperplanes (cutting planes) that are enforced as linear inequalities within the update subproblem. Because any subgradient defines the normal to a supporting hyperplane of the concave dual function, each such plane gives rise to the linear constraint:

It is added to the update problem. This prevents the update from leaving the validity range of the approximation. Constraints (Eq. 28) are referred to as bundle cuts and summarized as

In the remainder of the text. Yfantis et al. suggested computing the dual-variable update at every iteration by solving the following optimization problem [2, 3]: 

The constraint (Eq. 30b) imposes the identical trust-region bound employed in the BTM formulation (Eq. 17b). Enforcing all bundle cuts (Eq. 30c) from the outset can sometimes cause sluggish progress in the early stages of the algorithm. For this reason, the cuts are activated only when the current iterate lies within a user-specified distance from the anticipated optimum, controlled by a threshold. 

The resulting problem (Eq. 30) constitutes a quadratically constrained quadratic program (QCQP). Extensive numerical experiments indicate, however, that greater efficiency and robustness can be achieved by reformulating the bundle cuts using an epigraph representation: 

Here, the auxiliary variable 𝑣 serves as an upper bound on the current quadratic approximation (𝑡)(𝝀𝑁𝑝) at the candidate point (Eq. 32c). Thanks to this reformulation, the bundle cuts become simple linear inequalities (Eq. 32d), thereby improving numerical stability and accelerating overall convergence.

Preliminary numerical trials also indicated that the QNDA procedure generally reaches the region near the optimal point quite rapidly when applied to typical DMPC instances. Close to the solution, however, it frequently displays persistent oscillations that hinder satisfactory overall efficiency. To suppress this unwanted behavior, the mechanism for advancing the dual variables is changed—once the current iterate is sufficiently near the optimum—to a line-search procedure subject to constraints, rather than continuing with the trust-region approach:

The revised update rule (Eq. 33) essentially performs a step along the gradient direction computed from the quadratic model of the dual function. Constraint (Eq. 33b) guarantees adherence to the bundle cuts by adopting the epigraph-style representation from (Eq. 32c) to (Eq. 32d). The actual search direction appears in (Eq. 33c), and the allowable step length is restricted by (Eq. 33d), which takes the place of the previous trust-region bound. After determining this direction, the dual variables are incremented and then projected back into the nonnegative orthant to satisfy the nonnegativity condition imposed on them.

Activation of this constrained line-search update (Eq. 33) occurs only when the following condition is met:

With 𝜖𝑙 selected small enough to confirm that the iterates have indeed entered the immediate neighborhood of the optimal dual solution.

Figure 3 highlights the advantage offered by this line-search modification. The progression of primal and dual residual norms without line search is shown in Figure 3a. For comparison, Figure 3b shows the much smoother, faster decay of both residuals once the line-search update (Eq. 33) is engaged near the optimum.

Dual decomposition applied to dynamically coupled DMPC problems 

In the interest of completeness, this section briefly outlines the extension of dual decomposition to enable the distributed solution of DMPC problems with dynamic interconnections, as expressed in (4). To maintain brevity, coupling is restricted to the state equations alone, i.e., 𝐁𝑖𝑗 = 𝟎 for all 𝑖 ≠ 𝑗.

Begin by defining, for each subsystem 𝑖, the set of its direct neighbors N𝑖 ⊂ I, consisting of all indices 𝑗 for which 𝐀𝑖𝑗 ≠ 𝟎. Next, introduce auxiliary decision variables 𝐯𝑖𝑗 within subsystem 𝑖; these act as local replicas of the state trajectories from neighboring subsystem 𝑗. With this augmentation, the global centralized problem (4) can be recast in the equivalent form.

Under this rewritten structure (Eq. 35), the only remaining inter-subsystem links are the consistency (copy) constraints (Eq. 35f). Dual variables 𝝀𝑘𝑖𝑗 (𝑖 ≠ 𝑗) are then associated with each of these equality constraints. The resulting local Lagrangian for subsystem 𝑖 becomes

Thereafter, the complete MPC optimization task (Eq. 35) admits a fully distributed resolution: every subsystem independently minimizes its own Lagrangian. At the same time, the dual variables provide the necessary coordination across the network to enforce global consistency.

This section evaluates the performance of four distributed optimization algorithms—the subgradient method, the bundle trust method (BTM), ADMM, and the quasi-Newton dual ascent (QNDA)—by applying them to several carefully constructed DMPC benchmark instances. All test problems follow the general structure of problem (5), featuring a reference-tracking cost functional as in (2) supplemented by a terminal penalty term (3). 

The weighting matrices appearing in the objective were created randomly as symmetric positive definite matrices using the following procedure: 

Where the entries of 𝐍𝐱𝑖 and 𝐍𝐮𝑖 were independently sampled from the standard normal distribution, i.e., , 〖 ∈ (𝜇 = 0, 𝜎 = 1). 

System matrices 𝐀𝑖 were designed so that each resulting subsystem is asymptotically stable. For discrete-time linear systems, stability requires that every eigenvalue of 𝐀𝑖 lies strictly inside the unit disk [48]. Accordingly, each 𝐀𝑖 was formed via 

Here, 𝑖 is a diagonal matrix whose diagonal entries—the eigenvalues of 𝐀𝑖—were drawn from the continuous uniform distribution on (−1, 1), i.e., ∈(−1, 1). The elements of the orthogonal matrix 𝐒𝑖 were independently sampled from(−3, 3). The input-coupling matrices 𝐁𝑖 had entries drawn independently from the continuous uniform distribution  (−2, 2). 

State and input constraints for each subsystem were defined as the intersection of several ellipsoids centered at the origin: 

The shaping matrices 𝐆, and 𝐆_𝐮𝑖 were generated identically to the objective matrices (positive definite, random via the normal-based procedure above). The right-hand-side scalars 𝑝_𝐱𝑖, and 𝑝_𝐮𝑖, were independently drawn from  (1, 5). 

To ensure meaningful test cases, the desired reference trajectory 𝐱_𝑖^{ref,0∶𝑁𝑝} should itself be feasible. Its components were therefore sampled pointwise from  (−2, 2) at every time step until a feasible state (i.e., 𝐱_𝑖^ref, ∈ _𝑖. for all k) was obtained. The initial state 𝐱̃(𝑡₀) was generated following the same rejection-sampling procedure. 

Entries of the resource-interaction matrix 𝐑𝑖 were initially drawn from the uniform distribution  (1, 2). Subsequently, each entry was independently modified by multiplying it with a factor sampled from the discrete uniform distribution over {−1, 0, 1}: 

This random sign flip or zeroing guarantees that not every subsystem participates in every resource network and that some subsystems act as net producers. In contrast, others act as net consumers for any given resource. 

The upper bounds on resource usage were chosen deliberately so that the purely decentralized solution (corresponding to 𝝀 = 𝟎) becomes infeasible, whereas the fully coordinated (system-wide) problem remains feasible. To achieve this, the unconstrained centralized MPC problem (resource limits removed) was first solved to obtain the optimal input sequence. From this unconstrained optimum, the corresponding resource consumption profile was calculated: 

The resource bounds were then set by scaling down this unconstrained usage profile: 

After defining these tightened resource limits, the feasibility of both the decentralized problem (𝝀 = 𝟎) and the centralized (coupled) problem was explicitly verified to confirm the desired property.

The dimensions and counts of the subsystems were systematically varied across these ranges: 

- Total subsystems: 𝑁𝑠 ∈ {5, 10, 20, 50} 

- States per subsystem: 𝑛𝐱 ∈ {2, 3, 4, 5} 

- Inputs per subsystem: 𝑛𝐮 ∈ {2, 3, 4, 5}, subject to 𝑛𝐱 ≥ 𝑛𝐮 

- Shared resource channels: 𝑛𝐫 ∈ {2, 3, 4, 5}, subject to 𝑛𝐮 ≥ 𝑛𝐫 

- Length of prediction horizon: 𝑁𝑝 ∈ {10, 15, 20} 

All subsystems were given uniform internal complexity, so 𝑛𝐱𝑖 = 𝑛𝐱 and 𝑛𝐮𝑖 = 𝑛𝐮 held for every ∀𝑖 ∈ I. The quantity of local state constraints and input constraints per subsystem was deliberately matched to the number of states, giving 𝑛𝑐 = 𝑛𝐱. 

Ten separate benchmark instances were randomly created for every possible combination of the above parameters, producing an overall collection of 2400 different DMPC test cases. 

Both distributed algorithms and the individual DMPC subproblems were implemented in Julia [49], leveraging the JuMP algebraic modeling framework [50]. Every local DMPC optimization instance was processed by the commercial MIP solver Gurobi [51]. The trust-region update subproblem within BTM—linearly objective with affine equality constraints and convex quadratic inequality constraints—was also sent to Gurobi. Conversely, the quadratic update subproblems that arise in QNDA were solved using the nonlinear interior-point solver IPOPT [52]. 

Every simulation run was executed on a conventional laptop computer featuring an Intel(R) Core(TM) i5-6200U processor clocked at 2.30 GHz and equipped with 8 GB of main memory. 

Algorithm efficiency was quantified by estimating the total time required to complete each distributed optimization task according to the formula introduced in the study by Pflaum et al. [32]: 

Where 𝑁_iter stands for the total iteration count, 𝑇_comm is the assumed constant per-iteration communication latency between coordinator and subsystems, 𝑇_update^(𝑡)  is the coordinator’s CPU time for revising the dual variables in iteration 𝑡, and 𝑇sub^(𝑡),𝑖 denotes the solve time of local subproblem 𝑖 during iteration 𝑡. 

In a parallel distributed execution, the subproblems can run concurrently; however, the coordinator must wait for the last subsystem to finish before proceeding, so the slowest local solve bounds the effective per-iteration primal collection time. For consistency across all experiments reported below, the fixed communication delay was set to 𝑇comm = 800 ms. 

Parameter settings for DMPC problems 

The tuning parameters for the DMPC test cases were determined empirically through iterative experimentation to achieve an effective balance between numerical stability and convergence speed. 

For the subgradient method (SG), the initial step-size factor was fixed at (0) = 1; the same starting value was used for both the initial trust-region radius in BTM and QNDA. In all three cases, this parameter was subsequently adjusted following the rule. 

The ADMM scheme proved sensitive to overly aggressive tuning and frequently diverged under such conditions. Consequently, its starting penalty parameter was conservatively set to (0) = 10⁻³. The associated adaptation constants were chosen as 𝜏_incr = 1.25, 𝜏_decr = 1.1, and 𝜇 = 10. 

The threshold for activating bundle cuts in QNDA was placed at 𝜖𝑏 = 0.6 (Eq. 31), while the trigger for switching to the constrained line-search update (33) was set to 𝜖𝑙 = 10⁻² (see Eq. (34)). For both BTM and QNDA, the bundle retention depth (age parameter) was fixed at 𝜏 = 150. 

Every algorithm began with dual variables initialized to (0) = 𝟎; for ADMM, the auxiliary variables were likewise started at (0) = 𝟎 ∀𝑖 ∈ I. The approximate Hessian matrix required by QNDA was initialized as the negative identity matrix. 

Termination occurred after a maximum of 𝑡max = 500 iterations or whenever both the primal residual and the dual residual dropped below the common tolerance 𝜖𝑝 = 𝜖𝑑 = 10⁻². A complete overview of all selected parameter values appears in Table 1. 

Results for DMPC problems 

Performance outcomes for the distributed solution of the DMPC benchmark suite are visualized in Figure 4 and condensed in Table 2. Both the ADMM and QNDA algorithms successfully solved nearly all test instances. ADMM converged to a marginally higher proportion of problems. In contrast, QNDA demonstrated clear superiority in three key metrics: fewer iterations, shorter total computation time, and consistently lower residual values (better solution accuracy) among the converged cases. 

The classical subgradient method required a comparable number of iterations to QNDA for the problems it could solve, and its simpler update rule translated into lower per-iteration CPU cost, yielding better computational times. Nevertheless, it failed to converge to a substantially larger fraction of the benchmark set than ADMM and QNDA. The bundle trust method (BTM) delivered the weakest overall results across the DMPC test suite. 

Since the DMPC problems examined here are strictly convex, strong duality holds, and the optimal duality gap vanishes—the primal and dual objective values coincide at optimality. In contrast, nonconvex formulations (e.g., mixed-integer programs) typically exhibit a shrinking duality gap as the number of subsystems increases [53]. At the same time, a larger number of subsystems usually lessens the nonsmoothness of the dual function [3, 54]. 

For the convex benchmark problems considered in this study, these characteristics enable distributed solvers to deliver solutions that are noticeably more accurate overall, as shown in Figure 5. 

Importantly, the primal feasibility tolerance 𝜖𝑝 used to terminate the algorithms was intentionally kept fairly loose (Table 1). In most real-world control applications, extremely tight satisfaction of the coupling constraints is unnecessary. Moreover, techniques such as ADMM often converge to high-quality points quickly under moderate tolerances, whereas further tightening the tolerances demands a disproportionately large number of additional iterations [44]. 

Figure 5 shows that, across the converged runs, the average primal residual steadily decreases as more subsystems are included. Only successfully converged instances contribute to these averages. 

Figure 6 compares the average iteration count as a function of the number of subsystems for each algorithm. The subgradient method and BTM both exhibit a clear upward trend in the number of iterations required as the subsystem count increases. Conversely, ADMM, and especially QNDA, maintain low (and often decreasing) iteration counts, confirming their markedly better scalability and computational efficiency. 

Figure 7, on the other hand, indicates that varying the number of shared resources has comparatively little impact on the iteration counts of any algorithm. 

Figure 8 traces the convergence history of one representative test case involving three shared resources and a 15-step prediction horizon (yielding dual variables in R⁴²). Here, the subgradient method and QNDA reach convergence most rapidly. ADMM displays minor oscillations but still settles quickly. BTM, by contrast, suffers from large-amplitude oscillations and converges much more slowly relative to the other methods. The magnified view in Figure 8(b) reveals that these oscillations are concentrated near the optimal dual point, indicating that resource allocations are particularly sensitive to small adjustments in the dual variables. 

As reported in Table 2, the subgradient method achieves performance on par with QNDA across the cases it solves. This level of efficiency, however, depends on quite aggressive step-size choices that significantly reduce robustness. 

Figure 10 illustrates convergence curves for a more challenging benchmark with only two shared resources and a shorter 10-step horizon (𝝀 ∈ R¹⁸). Under these conditions, the subgradient method diverges wildly due to excessive oscillations triggered by the aggressive step size. All three algorithms converge successfully, with QNDA achieving the best speed. Notably, QNDA and BTM inherit the same trust-region radius as the subgradient method (and QNDA uses the same parameter to control its line-search steps). 

Figure 11 presents the terminal resource utilization for the problem in Figure 10. BTM, ADMM, and QNDA all achieve the optimal allocation. Figure 11(b) highlights that the subgradient method, however, stops at an infeasible point—over-allocating resource 2. 

This work has shown how dual decomposition techniques can be effectively employed to tackle distributed model predictive control (DMPC) problems in which subsystems are interconnected exclusively via shared scarce resources. The recently introduced quasi-Newton dual ascent (QNDA) algorithm was enhanced through two key modifications: an epigraph-based reformulation of the bundle cuts and the incorporation of a constrained line-search step applied once the iterates approach the optimal region. Extensive numerical experiments on a broad collection of benchmark instances confirmed the strong performance of this improved QNDA method. 

Although the alternating direction method of multipliers (ADMM) solved the largest proportion of test cases, the classical subgradient method delivered comparable efficiency to QNDA on the problems it solved, and QNDA consistently offered the most favorable trade-off between solution quality, iteration count, computational time, and numerical robustness. The benchmark results further revealed excellent scaling behavior for both ADMM and QNDA. As the number of subsystems increased, both the accuracy of the obtained solutions and the required number of iterations remained favorable. In contrast, the number of shared resources had only a marginal influence on iteration counts across all tested algorithms. 

Dual decomposition extends naturally to DMPC problems featuring dynamic coupling between subsystems. Future research could therefore focus on constructing a comprehensive suite of benchmark problems for dynamically coupled DMPC and on systematically comparing the performance of different distributed solvers in that setting. 

Linear dynamics governed all systems studied in the present paper. Because dual decomposition is an iterative, hierarchical procedure that requires multiple communication rounds before control inputs can be implemented, it is particularly well suited to processes characterized by slow dynamics—processes that can typically be approximated by linear models with acceptable accuracy over the prediction horizon. However, ongoing improvements in the speed and reliability of nonlinear optimization solvers are making nonlinear MPC increasingly attractive. Consequently, extending dual decomposition-based DMPC to systems described by nonlinear dynamics and nonlinear constraints represents a promising direction for further work. 

Nonlinear MPC problems are generally nonconvex, so convergence to an optimal dual solution does not automatically guarantee primal feasibility. In such cases, additional heuristic strategies are usually required to recover a feasible primal solution when applying dual decomposition to nonconvex subproblems. Nevertheless, even in the nonconvex regime, the dual function value remains a rigorous lower bound on the optimal primal objective. 

Finally, the literature still lacks realistic, industry-scale DMPC benchmark suites—particularly those involving a very large number of subsystems. While laboratory-scale examples are valuable for illustrating many theoretical and algorithmic properties, only application to truly large-scale industrial systems can fully reveal the practical benefits and limitations of a given coordination strategy.

None

None

None

None

Patrick O’Connor & Sean Murphy contributed to this work.

Authors and affiliations

Department of Power Systems Engineering, Faculty of Engineering, Trinity College Dublin, Dublin, Ireland
Patrick O’Connor & Sean Murphy

Corresponding author

Correspondence to Patrick O’Connor

Open Access The author(s) retain copyright. This article is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. It may be shared and adapted for non-commercial purposes with appropriate attribution, an indication of changes, and distribution of adaptations under the same license. Third-party material may be subject to separate terms identified in its credit line. View the license at https://creativecommons.org/licenses/by-nc-sa/4.0/.