A Modular Framework for Distributed Model Predictive Control of Nonlinear Continuous-Time Systems (GRAMPC-D)

The modular open-source framework GRAMPC-D for model predictive control of distributed systems is presented in this paper. The modular concept allows to solve optimal control problems (OCP) in a centralized and distributed fashion using the same problem description. It is tailored to computational efficiency with the focus on embedded hardware. The distributed solution is based on the Alternating Direction Method of Multipliers (ADMM) and uses the concept of neighbor approximation to enhance convergence speed. The presented framework can be accessed through Cpp and Python and also supports plug-and-play and data exchange between agents over a network.


Introduction
Model predictive control (MPC) is a modern control concept that attained increasing attention during the last decades [1,2] as it is capable to handle nonlinear systems while considering constraints on both states and controls. It is based on solving an optimal control problem (OCP) on a finite horizon and applying the first part of the control trajectory to the actual plant, corresponding to the sampling time ∆ t of the controller. At the next sampling instant, the horizon is shifted and the OCP is solved again. This iterative scheme is executed repetitively to stabilize the plant on an infinite horizon.
A main difficulty is the computational complexity of solving the OCP in real-time, which in turn requires an efficient implementation of suitable MPC algorithms. In the recent past, several toolboxes were published that provide adequate software frameworks such as ACADO [3] and ACADOS [4], VIATOC [5] or GRAMPC [6,7]. In case of distributed systems with a high number of controls and states, the classic centralized approach is not capable of solving the overall OCP in real-time anymore. Hence, algorithms for distributed model predictive control (DMPC) [8,9] have been in the focus over the last years. Their basic idea is to decouple the centralized OCP and to split it into multiple local OCPs that can be solved in parallel. The expectation is to compensate the higher computational complexity due to the decoupled formulation as well as the increased communication effort by the parallel structure. There are multiple approaches to distributed algorithms for optimal control problems, such as sensitivity-based algorithms [10], the augmented Lagrangian based alternating direction inexact Newton method (ALADIN) [11,12] or the alternating direction method of multipliers (ADMM) [13] that is also used in the presented framework.
The difficulty of an efficient implementation is drastically higher in case of DMPC than for classic MPC algorithms, as a potentially high number of subsystems, so-called agents, have to be managed. Several toolboxes for DMPC have been published as well. Linear discrete-time systems are considered in the DMPC-Toolbox [14] that is implemented in Matlab. The PnPMPC-TOOLBOX [15] focuses on the plug-and-play functionality and provides an implementation in Matlab that considers continuous-time and discrete-time linear systems. Several algorithms are implemented in the Python-Toolbox DISROPT [16] regarding distributed optimization problems. ALADIN-α [17] is the most recent published toolbox that provides a Matlab implementation of the ALADIN algorithm. However, there is a lack of a DMPC framework that provides an implementation tailored to embedded hardware with the focus on real-time capable distributed model predictive control. Many real-world problems such as smart grids or cooperative robot applications are only equipped with weak hardware on the subsystem level that is not able to handle complex computation tasks in an appropriate time. Hence, for realizing distributed controllers on actual plants, an implementation optimized on execution time is required to enable real-time control. Furthermore, providing the possibility of communication between agents over a network is essential for a DMPC framework designed to control actual plants. The restriction to neighbor-toneighbor-communication decouples the agents communication effort from the overall system size by bounding it to the cardinality of its neighborhood. The focus on real-world plants requires the system class to cover nonlinear dynamics including couplings between the agents in both dynamics and constraints.
The presented framework, in the following GRAMPC-D, provides an opensource C ++ implementation that is capable of solving optimal control problems in a distributed manner with a per-agent computation-time in the millisecond range. The underlying minimization problems are solved with the MPC toolbox GRAMPC that is suitable for embedded hardware implementations. However, other toolboxes for solving the local minimization problem can be used as well. To enable actual distributed optimization, a socket-based TCP communication is provided to allow agents to exchange data over a network. Furthermore, a Python interface is provided in addition to the C ++ -interface using the software module Pybind11 [18]. The Python interface combines both the functionality of Python and the performance of C ++ as it only wraps the C ++ interface while the actual code executing the DMPC algorithm is still running in C ++ . Furthermore, it allows for fast and efficient prototyping when developing a controller for a distributed system as both a centralized as well as a distributed controller can be derived based on the same problem description. The convergence behavior of distributed controllers can be improved by optionally using the concept of neighbor approximation. Thereby, the generated problem description of each agent is adapted to additionally approximate parts of its neighbors OCP and by this to improve the solution of its local OCP in each iteration, leading to an enhanced convergence behavior of the overall algorithm. The modular structure of GRAMPC-D enables modifying the overall system in the sense of plug-and-play by including or removing agents or couplings at run-time. Supporting plug-and-play features is a core functionality for a DMPC framework with focus on embedded systems, as the assumption of a static system description does not hold for a large number of real-world plants.
The paper is structured as follows. Section 2 outlines the considered class of coupled systems and OCP formulation. The DMPC framework GRAMPC-D is introduced in Section 3 including the ADMM algorithm as the method of choice for the algorithm. In addition, the concept of neighbor approximation is explained and the implemented algorithm for the crucial task of penalty parameter adaption is presented. The modular structure of the framework is presented in Section 4. Finally, simulation examples in Section 5 show the effectiveness and modularity of the DMPC framework, before conclusions are drawn in Section 6.
Throughout the paper, each vector x ∈ R n is written in bold style. Standard p-norms 1 p will be used as well as the weighted squared norm defined by x 2 P = x T P x with a positive (semi-)definite square matrix P . The stacking of individual vectors x i , i ∈ V from a set V is denoted by x = x i i∈V . As far as time trajectories are concerned, the explicit dependency on time t may be omitted to ease readability. The derivative with respect to time is written using the dot notationẋ(t) = d dt x(t).

Problem Description
The presented DMPC framework considers multi-agent systems that can be described by a graph G = (V, E) with the sets of edges E and vertices V. Each vertex represents an agent, while each edge between two vertices stands for a coupling between the corresponding agents. The couplings may be both uniand bi-directional. The considered optimal control problem for the coupled nonlinear system is given by states x i (t) ∈ R nx,i , controls u i (t) ∈ R nu,i and the horizon length T > 0. Each agent may have a general nonlinear cost function l i : R nx,i × R nu,i × R → R and terminal cost V i : R nx,i × R → R. The overall cost function (1a) is given by the sum over the individual cost functions. The subsystem dynamics (1b) are defined by the functions f i : R nx,i × R nu,i × R → R nx,i and f ij : The OCP additionally considers nonlinear equality constraints (1d)-(1e) and inequality constraints (1f)-(1g) with the functions g i : ij as well as box constraints (1h) for the control input u i of each agent i ∈ V. The neighborhood N i of agent i ∈ V is given by two sets that differ in the direction of the coupling, sending neighbors N ← i and receiving neighbors N → i , see Figure 1 for an example. States and controls of sending neighbors have an explicit influence on the dynamics of the agent i in form of functions f ij , see (1b). Receiving neighbors are neighbors of agent i that are explicitly influenced by this agent, hence states or controls of the agent are part of a function f ij of receiving neighbors. While a neighbor can be both, receiving and sending, this separation is going to be beneficial in the ADMM algorithm by reducing unnecessary computation and communication effort.
The dynamics (1b) of each agent i ∈ V are neighbor affine in the sense that the dynamics consists of a function f i that depend only on states and controls of the agent and a sum of functions f ij that depend on states and controls of the agent and one neighbor. The constraints (1d)-(1h) on each agent are separated into constraints that depend on states and controls of the agent, given by g i , h i and the box constraints (1h), and constraints g ij and h ij depending on states and controls of the agent and one neighbor, similar to the dynamics.
The considered OCP formulation (1) covers a wide class of distributed systems, e.g. cooperative transport [19] and scalable systems such as smart grids [20]. This generic system description combined with the focus on a timeefficient implementation opens a wide spectrum of usability for the presented DMPC framework.

Distributed model predictive control
Optimal control problems for coupled systems as in (1) contain a large number of states and controls. This leads to a significant computational effort that is challenging for standard MPC algorithms to be handled in real-time. DMPC algorithms instead assume that each of the distributed subsystems are equipped with a dedicated control unit that is capable of solving a reduced optimal control problem. The idea based on this assumption is to decouple the global OCP and spread the computation effort over the set of agents in parallel. Overlying algorithms ensure convergence of the local solutions to an optimal solution for the overall system. While the computational complexity and communication effort of algorithms for DMPC is higher than solving the central problem, the expectation is to compensate this disadvantage by the parallel structure. In the presented DMPC framework, the well-known ADMM algorithm [13] is employed in a continuous-time setting [21].

ADMM algorithm
The ADMM algorithm enables to spread the computation effort of the global OCP (1) completely on distributed agents. As a starting point, the global OCP (1) is brought into a decoupled form for each agent i ∈ V by introducing local copiesx ji (t) ∈ R nx,j andū ji (t) ∈ R nu,j for the states x j and controls u j of each sending neighbor j ∈ N ← i , i.e. min w,z i∈V These local copies (x ji ,ū ji ) represent new control inputs for the agent i and can be seen as a proposal of agent i for its neighbors j ∈ N ← i . Equivalence between the local copies and the original variables is ensured by introducing the consistency constraints (3i) and (3j) with the coupling variables z x,i (t) ∈ R nx,i and z u,i (t) ∈ R nu,i . In (3) and the following, the notation is used. The ADMM method is based on the Augmented Lagrangian formulation [21,22]. Regarding the continuous-time setting used in this paper, the consistency constraints (3i) and (3j) are accounted for in the cost functional To ease notations, the multipliers and penalty parameters are stacked according to The corresponding dual problem to (3) can be written as with the primal variables (w, z) and the dual variables µ. The ADMM algorithm solves the max-min-problem (7) by repetitively executing the three steps with the iteration counter q. The minimization with respect to the coupling variables z (8b) can be solved analytically while the steepest ascent is used in (8c). Important to note is that each step can be subdivided into fully decoupled steps for either agent i ∈ V. Hence, the algorithm is fully distributable which allows to spread the computation effort over all agents. The resulting ADMM algorithm for each agent is given in Algorithm 1. It consists of the computation steps 1, 3, 5, the communication steps 2, 4, 6, and the evaluation of a convergence criterion in Step 7. The algorithm starts with

Algorithm 1 Alternating direction method of multipliers
2: Send local copiesx q ji andū q ji to sending neighbors j ∈ N ← i 3: Compute coupling variables 4: Send coupling variables z q i to receiving neighbors j ∈ N → i 5: Compute Lagrange multipliers 6: Send Lagrange multipliers µ q i to sending neighbors j ∈ N ← 8: STOP 9: else 10: set q = q + 1 and go to Step 1. 11: end if an initialization of corresponding variables. The local OCP (9) is minimized in Step 1 with respect to the local variables w i . This minimization represents the main computation effort of the overall algorithm. In Step 2, the trajectories of the local variables are sent to the sending neighbors j ∈ N ← i of each agent i ∈ V. The analytic solution for the minimization with respect to the coupling variables (8b) is given in Step 3 of the ADMM algorithm, before they are sent to the receiving neighbors j ∈ N → i of each agent in Step 4. The third computation step is given in Step 5 by a maximization with respect to the Lagrange multipliers µ. In Step 6, the result of the maximization step is sent to the sending neighbors j ∈ N ← i of each agent i ∈ V. A convergence criterion is checked in Step 7. If it is satisfied or the iteration counter has reached its maximum, the algorithm stops and returns the current trajectories. Otherwise, the iteration counter is increased and the algorithm returns to Step 1.

Neighbor approximation
In practice, the convergence speed of the ADMM algorithm can be enhanced by anticipating the actions of the neighbors in the own agents optimization. The concept of neighbor approximation was introduced in [19] and extended in [23] and relies on the neighbor affine structure of the dynamics (1b)-(1c) and constraints (1d)-(1h). The basic idea is to use the already introduced local copiesx ji andū ji to approximate parts of the neighbors OCP. The expectation is that the additional information about the neighborhood improves the local solution of each agent and thus the convergence behavior of the overall algorithm. This also reduces the number of required ADMM iterations until convergence is reached, which has been confirmed in numerical evaluations in [19] and [23]. In practical experience, the reduced number of ADMM iterations can compensate for the increased complexity of the extended OCP which can lead to a significantly decreased computational effort [23].
The neighbor approximation implemented in GRAMPC-D is modular in the sense that the neighbors cost, constraints, dynamics and each combination of the three can be considered.

Neighbor cost
The global cost to be mininized (2) consists of the single cost functions of the agents i ∈ V. The local copies of the neighbor variablesx ji andū ji , j ∈ N i , can be used to anticipate the neighbors cost J j (ū ji ; x j,0 ) on the local level of agent i ∈ V, i.e.
The normalization with the factors is necessary in order to avoid that the neighbors cost function would appear in the overall cost function multiple times. Approximating the neighbor costs is especially beneficial in examples with a strong dependency on the other agents costs and enables the agent to anticipate the neighbors control action to minimize its local costs. It is recommended to combine the neighbor cost approximation with the approximation of the neighbor dynamics introduced in the following lines.

Neighbor dynamics
Similar to the neighbor cost consideration, the neighbor affine structure of the dynamics (1b)-(1c) can be exploited to approximate the neighbor dynamics and therefore to improve the quality of the local copiesx ji andū ji . To this end, the local dynamics (1b)-(1c) are extended by the approximate neighbor dynamicṡ with the initial conditionx ji (0) = x j,0 . The dependencies of the neighbor's states and controls in f j and f ji are decoupled using the local copiesx ji andū ji . However, it is not possible to decouple further functions f js as these depend on states and controls of agents s for which agent i has in general no local copies. For consistency, the external influence is introduced with v ij (t) ∈ R nx,i that captures the remaining terms of the neighbors dynamics. The external influence is considered in the approximated neighbor dynamics (14) by introducing local copiesv ji (t) ∈ R nx,j . Thereby, the whole dynamics of neighbor j is approximated in (14). To ensure convergence of the local copiesv ji to the original variables v ij , the consistency constraints are introduced and replace the consistency constraints in (3i)-(3j) regarding the states. Note that the local copies of the statesx ji are not considered as control variables anymore, but are determined by the differential equation (14). Instead, the local copies of the external influencev ji serve as new local control variables. In summary, the stacked notations (4) and (6) are adapted according to and w = w i i∈V , z = z i i∈V , µ = µ i i∈V (18) with Lagrangian multipliers µ v,ij (t) ∈ R nx,i , µ v,ji (t) ∈ R nx,j , coupling variables z v,ij (t) ∈ R nx,i , and penalty parameters ρ v,ij (t) ∈ R nx,i , ρ v,ji (t) ∈ R nx,j .

Neighbor constraints
In addition to the consideration of the neighbor cost and dynamics within the local OCP of agent i ∈ V, the constraints (1d)-(1h) of each neighbor j ∈ N i of agent i ∈ V can be taken into account by adding to the local OCP (9). Again, the constraints are decoupled from the neighbors states and controls x j and u i by using the local copiesx ji andū ji . As discussed before, this concept is restricted to the agent constraints of each neighbor j ∈ N i and the coupling constraints between neighbors j and agent i, while further coupling constraints between neighbor j and its neighbors s ∈ N ← j \ {i} depend on states and controls of agents s for which in general agent i has no local copies.

Penalty parameter adaption
The update of the penalty parameters in the matrices C i and C ji in (7) is crucial for a fast convergence of the ADMM algorithm. The adaptation method implemented in GRAMPC-D follows a proposal in [13, Section 3.4.1] for the optimization problem The proposed adaption algorithm is given by with the primal residual r q = Ax q +Bz q −c, the dual residual s q = ρA T B(z q − z q−1 ). The basic idea is to keep both within a factor of µ of one another. Following this idea for the OCP (7), the primal and dual residuals are given by

3:
if γ q > γmax then 4: γ q = γmax 5: else if γ q < γ min then 6: γ q = γ min 7: end if 8: else 9: γ q = 1 10: end if 11: ρ q = γ q ρ q−1 To reduce the number of tuning parameters, µ = 1 is chosen which results in an equality instead of the inequality in (21). To further simplify the implementation, the equality is evaluated element-wise and at each discrete time step δ k = T N −1 with N as discretization of the predicted horizon. This results in the condition for each discrete time step δ k and the norm evaluated element-wise. The condition (23) can be reformulated in form of the update law with m as index for an arbitrary element in (23). The implementation is presented in Algorithm 2. At first, the division through small numbers, especially zero, is caught to prevent numerical issues. The factor γ q is computed afterwards and bound between γ min and γ max , before the new penalty parameter is calculated by ρ q = γ q ρ q−1 .

Modular framework
GRAMPC-D is implemented in a modular fashion in order to achieve a scalable and flexible implementation. At first, the modular structure is explained before the capability for plug-and-play scenarios is laid out.

Modular structure
The main parts of GRAMPC-D and their interaction are presented in the following. Due to the modular concept, the implementation of GRAMPC-D can be subdivided into single modules that are composed depending on the chosen type of controller, centralized or distributed, as the structure of GRAMPC-D differs between the two cases. Both are visualized in Figure 2. Either structure is generated automatically by choosing the corresponding controller type without further required interaction of the user. Both the distributed and centralized structure are scalable due to the modular concept and therefore suitable to handle large or complex systems. The distributed control structure in the left part of Figure 2 assumes that each agent only has access to its local variables and communication is required to acquire data from other agents. Thus, the central part in the distributed setup is the communication interface. While it enables to exchange data between agents, the actual implementation depends on the chosen type of communication interface. If the DMPC is simulated on a single processor, there is no need to actually send data over a network. Instead, a central communication interface is provided that exchanges data pointers, which is a significant difference in performance. If the ADMM algorithm is implemented in an actual distributed setup, each agent creates its own local communication interface that enables to exchange data over a network. The corresponding protocol is encapsulated into the local communication interface due to the modular concept that enables implementing multiple protocols and switching between them. In either case, each agent creates a local solver that contains the local OCP depending on the neighborhood and the chosen optimization parameters such as neighbor approximation. Hence, tasks like decoupling the global OCP by introducing local copies are done automatically in the background. The ADMM algorithm is implemented inside the local solver with an abstract implementation of the minimization problem with respect to the local variables. This enables implementing multiple solvers and switching between them without changing other parts of the software, although GRAMPC is chosen as default. The remaining two important modules are the coordinator and the simulator. The ADMM algorithm assumes a fully synchronized execution, which has to be guaranteed even in a distributed setup. This synchronization is handled by the coordinator by triggering each step of the algorithm and waiting for a response of each agent before sending the following trigger. The last module is an integrated simulator that enables simulations independent of the chosen controller or the specific system.
A more simple structure is generated if a centralized controller is chosen, see right part of Figure 2. In this case, the global OCP (1) is solved in a centralized manner including all agents dynamics and having knowledge of all variables. The centralized setup implicitly synchronizes the execution of the algorithm without the need for a coordinator. The only remaining module is the simulator that is used to simulate the overall system.
Each part of GRAMPC-D is interchangeable by alternative software. As already mentioned, the MPC toolbox GRAMPC is used by default to solve both the global OCP (1) in case of a centralized controller and the underlying reduced OCP (9) in Step 1 of the ADMM algorithm in case of a distributed controller. GRAMPC is tailored to embedded hardware and by this a natural choice, but if another solver is desired, e.g. with a stronger focus on precision instead of computational speed, then the only required change is to overload the class regarding the solver with a new implementation. The same holds for each part of the framework such as the implemented communication protocol. The TCP protocol is provided by default, but alternative protocols can be implemented by overloading the local communication interface.

Rapid prototyping
Both the local and global OCPs are dynamically generated at run-time based on the same problem description provided by the user, see Figure 3. In case of a centralized controller, the global OCP is generated by the central solver while in case of a distributed controller, the local OCP is generated for each agent individually based on its neighbors and the optimization parameters. Therein, the corresponding flags for neighbor approximation are defined that state whether additional variables have to be initialized, e.g. local copiesx ji andū ji or the external influence v ij , see Section 3.2. This results in a convenient prototyping process while designing controllers, as each type of controller can be automatically generated and evaluated based on the same problem description. To further support the efficiency of prototyping, the possibility of multi-threading can be activated to spread the computation effort on each available core of the processor and thus to speed-up the computations.

Plug-and-play functionality
Plug-and-play is a core feature for the usability of a framework in the field of DMPC, see e.g. [24,25,26], where the OCP structure and size may change dynamically due to the removal or plug-in of agents. The generation of the OCPs (9) during run-time and the modular structure of the framework allows to integrate plug-and-play functionality in the network. If an agent enters the system in the distributed control setting (left part of Figure 2), the coordinator informs the direct neighbors of the agents to include the corresponding variables. Hence, only the OCPs (9) of the direct neighbors are updated and the new agent is integrated into the network. If an agent leaves the network, the  agents and the coordinator delete the agent from their internal lists of active agents. This results in a smooth transition between two problem formulations.

Simulation examples
The modular framework GRAMPC-D is evaluated for different examples. The scalability is shown for a coupled spring-mass system. The plug-and-play functionality and the concept of neighbor approximation are demonstrated for a smart grid and a coupled water tank network, respectively. Finally, a distributed hardware implementation for a system of coupled Van der Pol oscillators is considered by communicating over an actual network using the TCP protocol.
In each example, the terminal and integral cost in (2) are chosen quadratically with the positive (semi-)definite weighting matrices P i ∈ R nx,i×nx,i , Q i ∈ R nx,i×nx,i and R i ∈ R nu,i×nu,i . The desired state to be controlled is given by x i,des . The computation times are measured on an Intel i5 CPU with 3.4 GHz using Windows 10. The communication effort is neglected if not stated otherwise.

Scalable system
The scalability of GRAMPC-D is shown for a system consisting of a set of masses that are coupled by springs. Each mass is represented by an agent i ∈ V and is described by the differential equations with the position (p x,i , p y,i ) of the respective mass in the x− and y-axis and the respective controls (u x,i , u y,i ). This results in the state and control vectors The spring is relaxed at the length δ 0 = 1 m. The spring constant is given by c = 0.5 N m −1 and each agent has a mass of m i = 7.5 kg. The function δ ij (p x,i , p y,i ) = (p x,i − p x,j ) 2 + (p y,i − p y,j ) 2 computes the distance between two agents i and j. The dynamics (26) can be split into functions f i and f ij corresponding to the neighbor-affine form (1b). The weighting matrices are set to (SI units are omitted for simplicity) [5,2,5,2] , R i = diag [0.01, 0.01] (28) with the desired state It was shown in [27] that the computation time per agent is nearly independent of the system size, whereas in the central MPC case the computation time rises drastically. The simulation results for a system with 40 × 40 agents are given in Figure 4. It can be seen that the trajectories of the cost are quite similar. While the distributed solution is slightly suboptimal, it would converge to the centralized solution by increasing the number of ADMM iterations. The computation time for each time step, however, is 1072.58 ms for the centralized controller, while the distributed controller requires a maximum of 9.59 ms and an average of 2.23 ms per agent.

Plug-and-play
The plug-and-play capability of GRAMPC-D is presented using an exemplary setup of a smart grid. The network is described by a set of coupled agents that represent non-controllable power sinks and sources, such as private households, industry or renewable energy as well as controllable power plants. The dynamical behavior of the agents is generalized by describing them as generators with a mechanical phase angle θ i (t) ∈ R that may differ from the phase of the grid. The corresponding dynamics with friction constant κ > 0 and the moment of inertia I is given in a Power plant Power sink Power sink Fig. 5 The plug-and-play capability of GRAMPC-D is presented using the simulation of a smart grid. At the beginning of the simulation, only one power sink is connected to the power plant. During run-time, the second power sink is connected to the first one, i.e. the power plant has to supply both using the same connection.
neighbor-affine form and describes the dynamical behavior of the phase shift φ i (t) ∈ R given by between the phase Ωτ of the grid with frequency Ω and the mechanical angle θ i , see [28]. Hence, the state vector is given by In (30), P source,i describes the generalized non-controllable power, e.g. the demanded power for private households and industry or the generated power by renewable energy. The coupling between two agents consists of the maximum transferable power P max,ij that depends on the phase shift angle φ j − φ i between agent i ∈ V and its neighbors j ∈ N ← i . Agents that describe power plants have a controllable input u i that is used to stabilize the grid. A normalized parameterization is used with Ω = 1 Hz and I = 1 J s 2 , the friction term set to κ = 1 × 10 −3 J s and the maximum transferable power to P max,ij = 0.1 J s −1 . The weighting matrices are set to with the desired state The first element of the desired state x i,des is set arbitrary, as there is no desired value for the phase and the first state is not weighted in the cost functional.
The implemented setup is visualized in Figure 5. At the start of the simulation, one power plant is given that supplies one non-controllable power-sink such as a private household. During run-time, an additional power sink is coupled to the first one, i.e. the power plant has to supply both using the same connection. The simulation results are shown in Figure 6, starting with the first power sink that is connected to the power plant. It can be seen that the phase difference between the power plant and the household converges to a stationary value that leads to a transmission of the demanded power. Furthermore, the angular velocity of the phase shift converges to zero. At simulation time t = 20 s, the second power sink is plugged in, leading to an additional power demand. Consequently, the phase shift between the power plant and the first power sink increases and the additional power is transmitted. The phase shift between the first and second power sink is adapted accordingly. The computation time for the distributed controller is given by a maximum of 84.67 ms and an average of 4.78 ms per agent using a step size of ∆ t = 100 ms. The average time is significantly lower due to the convergence criterion of the ADMM algorithm.
This plug-in and plug-out functionality of agents is supported at any moment during the simulation even if the controllers run on distributed hardware and communicate over a network. GRAMPC-D can handle planned changes in the system such as shown in this simulation example as well as spontaneous disconnections due to a broken network connection.

Neighbor approximation
The concept of neighbor approximation is evaluated for a system of water tanks that are coupled by pipes, see Figure 7. Only the first water tank has a controllable input u 1 (t) ∈ R. The last water tank has a constant outflow d 5 = 0.01 m 3 s −1 and the desired water height x 5,des = 3 m. In addition, the inequality constraint h i ≤ 3 m of a maximum water height has to be satisfied by all tanks. The dynamics of each water tank is given bẏ with the water height h i (t) ∈ R and the state vector x i = h i . The area of each water tank is given by A i = 0.1 m 2 and the diameter of the pipes by a ij = 0.005 m 2 . The weights for the cost functions are set to The simulation is run with a distributed controller both with and without neighbor approximation using the same set of parameters for the ADMM algorithm and GRAMPC. The convergence of the ADMM algorithm is shown in Figure 8 for both simulations. The simulation with neighbor approximation converges smoothly to the optimal solution and satisfies the convergence criterion after 7 iterations while 89 ADMM iterations are required without neighbor approximation. Note that the cost is rising instead of falling as the solution is infeasible until the algorithm converges. The improved convergence behavior with neighbor approximation comes with a higher computational complexity per ADMM iteration. However, if a convergence criterion is used, the decreased number of ADMM iterations per time step compensate for the higher computational complexity. The required computation time for the 89 ADMM iterations without neighbor approximation is 955.6 ms and for the 7 iterations using neighbor approximation 203.8 ms. Note that the same configuration for the ADMM algorithm and GRAMPC is used in this evaluation to provide a comparable result. The standard ADMM algorithm may converge within less iterations if more computation effort is spent per iteration while the 1 2 3 4 5 Fig. 7 The concept of neighbor approximation is shown at a simulation example of coupled water tanks. Only the first one has an input and only the last one has a desired water height while being disturbed by a constant outflow.  algorithm may require even less time if the parameters are tuned for neighbor approximation.

Distributed hardware implementation
This section shows the capability of GRAMPC-D to solve the ADMM algorithm on distributed embedded hardware. The following simulation is run on four Raspberry Pi 3B+ using Raspberry Pi OS that are connected via Ethernet. Each agent runs on an individual Raspberry Pi while the coordinator and the simulator use the same hardware. The local communication interface from GRAMPC-D based on the TCP protocol is used. The simulation example consists of three coupled Van der Pol oscillators [29] with state p i (t) ∈ R, control u i (t) ∈ R, the uncoupled oscillator constant α 1 = 1 m −1 s −2 , and the coupling constant α 2 = 1 s −2 . The state vector and desired state are given by x i,des = 0 m 0 m s −1 T (38b) with the weighting matrices set to The simulation is run for 10.000 time steps using a fixed number of q max = 5 ADMM iterations. Figure 9 shows the required time for computation and communication in each time step. The computation time to execute 5 ADMM iterations amounts to 10 ms per agent while the average time required to solve the ADMM algorithm in a distributed manner using the TCP protocol for the communication is 66.51 ms. Hence, the average effort for the communication is 56.51 ms including 72 communication steps. Since all agents have to be synchronized for the ADMM algorithm, either of them has to wait for the slowest agent at each step of the algorithm. All five ADMM iterations can be executed in 48 ms with the worst case time of 114 ms. This results in the minimum time for each communication step of 0.53 ms, an average time of 0.79 ms and a maximum time of 1.44 ms. This time includes preparing the data to be sent, sending and receiving it and recreating the sent data structure from the byte array. These are plausible values as a ping to the loopback address 127.0.0.1 already requires 0.14 ms in average and 0.22 ms at maximum. These results show that the main effort in distributed optimization is the communication effort that requires 82.3% of the overall time in average for this simulation example.

Conclusions
The open-source, modular DMPC framework GRAMPC-D is presented in this paper that enables to solve scalable optimal control problems in a convenient way and to stabilize plants using distributed model predictive control. This problem description can be used for both a centralized and a distributed controller. In the distributed setting, the global optimal control problem is automatically decoupled and solved in a distributed manner using the ADMM algorithm. The convergence behavior of the ADMM algorithm can be improved by sugint he concept of neighbor approximation that allows the agents to anticipate the actions of their neighbors. The presented DMPC framework supports plug-and-play to connect and remove agents at run-time. Besides solving the ADMM algorithm on a single processor, it is possible to solve the local optimal control problems on distributed hardware. The local communication interface enables communication between agents over a network using the TCP protocol. By default, GRAMPC-D uses the MPC toolbox GRAMPC for solving the local optimal control problems on agent level, which is suitable for real-time and embedded implementations. GRAMPC-D is licensed under the Berkeley Software Distribution 3-clause version (BSD-3) license. The complete source-code is available at Github https://github.com/grampc-d/grampc-d. Future work will use the modular structure to extend GRAMPC-D. For example, communication protocols besides TCP can be provided or alternative solvers to GRAMPC for the underlying minimization problem implemented to increase the usability and flexibility of the framework.