Agent-Based Decentralized Model Predictive Control for Plants With Multiple Identical Actuators

This article proposes a decentralized model predictive control (DMPC) algorithm without communication for systems consisting of multiple identical, independent actuators acting on a single central plant. The particular system design is relevant for applications where modularity is paramount and for highly dynamic systems, e.g., battery emulator systems controlled by multiple dc–dc converter modules. Each agent, consisting of an actuator and a DMPC, controls a virtually scaled version of the plant to implicitly consider the effects of other agents. The set of DMPCs achieves the same plant performance as a corresponding centralized model predictive controller (CMPC) in unconstrained operation. Also, the states of the independent agents converge toward the globally optimal CMPC solution. This is obtained by dividing the CMPC’s objective function into local objective functions related to the subsystems. The optimality and stability of the DMPC in unconstrained operation are shown analytically. The stability of control input-constrained operation is analyzed by computing the region of attraction. Numerical studies of a battery emulator system compare the performance of the DMPC with the global optimum in detail.


I. INTRODUCTION
D ECENTRALIZED model predictive control (DMPC) has been investigated intensively since the early 2000s [1], [2], [3]. It offers the possibility to apply model predictive control with its strengths to systems that cannot be controlled by a centralized instance. Two of the main issues of centralized model predictive control (CMPC) are: 1) the Manuscript  complexity of solving centralized problems grows with the number of considered states and 2) measurement data of the entire system must be communicated to the central processing unit and control signals must be transmitted back to the decentralized actuators [4]. Both aspects become even more critical if highly dynamic systems involving short sampling intervals are treated. To overcome these obstacles, DMPC approaches divide the centralized control problem into several less complex tasks, each of which is solved by a local controller. In general, the cost of this decentralization is the loss of performance with respect to the CMPC. Therefore, most of the DMPC approaches introduce different kinds of communication between local controllers in order to improve the global performance [5]. However, noncooperative communication strategies, where each local controller optimizes a local objective function, do not necessarily improve the control performance and might even cause closed-loop instability [6], [7]. Moreover, communication requires additional infrastructure and can also be technologically restricted, e.g., in earlier mentioned highly dynamic systems [8].
Decentralization is based on decomposing the overall system into smaller subsystems. The majority of publications on DMPC focus on networked systems comprised of interconnected subsystems [9], where the decomposition of the overall system is rather clear. However, in some applications, there are systems with fundamentally different structure, where a single plant or central system is controlled by several self-contained, dynamic actuators. Due to the completely different roles of the actuators and the central system in the overall system, the formulation of subsystems is not straightforward.
Addressing a specific class of those applications, the main contribution of this work is a DMPC algorithm without communication that achieves global CMPC optimality regarding central system states in unconstrained operations. Each actuator is equipped with a DMPC, which is together defined as an agent (see Fig. 1). The proposed algorithm targets the aforementioned applications, where a central system is controlled by multiple actuators. All actuators are identical, work independently of each other, and are not directly coupled. However, states of the central system generally imply a feedback effect on actuator states, which means that the actuators are indirectly coupled through the states of the central system (see Fig. 1). This makes controlling the overall system in a decentralized manner challenging.
The special system design with identical actuators is often found in the literature as it involves some important advantages, such as improved performance, modularity, scalability, efficiency, and cost. Compared to a single actuator, multiple smaller actuators are potentially able to respond faster to rapidly changing loads as parallel-connected dc-dc converters show [10]. Modular design allows to adapt easily to changing plant requirements by just adding or removing actuators. Under this perspective, a control strategy for parallel-connected inverters for uninterruptible power supplies is presented in [11]. Also, efficiency during operation benefits from modularity. Adjusting the number of active actuators depending on the actual load is generally better than running a large actuator under partial load. This can be seen in [12], where an optimal control strategy for redundant refrigeration circuits acting on a cooling chamber is presented. Similarly, multiple parallel compressors can be installed in a compressed air system [13]. From an economic point of view, using multiple identical actuators "off the shelf" is more efficient than designing a single one that fits the scale of the application [14].
The overall control problem of the investigated class of systems comprises two objectives, namely, load sharing among the participating actuators and controlling the central system output. In particular, control methods for parallelconnected dc-dc converter systems are widely represented in the literature. Besides the conventional droop control, various active current sharing approaches have been proposed, including sophisticated control concepts [15], [16], [17]. DMPC approaches without communication between local controllers have also been reported: The method presented in [18] aims at parallel converters controlling the voltage at a central capacitor and particularly focuses on converter interleaving, a parallelconverter-specific topic. In [19], subsystems for the DMPC are composed of the central system, the local actuator, and a fictitious actuator representing all other actuators. In contrast, each agent controls the voltage at its own local capacitor in [20], i.e., no actual central system exists, and the subsystems are therefore directly derived.
An alternative decentralized control framework often found in practice and literature is the decentralized proportionalintegral-derivative (PID) controller. Although numerous methods for designing decentralized PID controllers for multivariable systems have been developed, the implementation is rather complicated, requires high tuning effort, or is simply based on heuristics [21], [22], [23]. Moreover, loop interactions can cause sacrifices in performance or even lead to instability and therefore impede the implementation of this control concept [24], [25]. Another drawback of decentralized PID controllers is that predictive consideration of constraints is not directly possible.
This work proposes a DMPC method that is based on creating virtual subsystems by assigning a scaled version of the central system to each agent to implicitly consider the effects of other agents. The method works without communication, and to set up local controllers, the number of agents is the only necessary information. The local objective functions are based on an exact division of the CMPC objective function. This allows to compare the performances of the proposed DMPC and the corresponding, equally weighted CMPC in a fair manner. Although agents do not communicate, the proposed DMPC performs with the global optimality of the CMPC regarding all central system states in unconstrained operation, and the agent states converge to the solution of the CMPC. In control input-constrained operation, the DMPC deviates from the CMPC optimality but stays close to it and converges when control input constraints become inactive.
The remainder of this article is structured as follows. In Section II the system description is presented. A detailed description of the DMPC is given in Section III, and in Section IV, the structure of the CMPC can be found. The optimality and stability of the approach are analyzed in Section V. Section VI demonstrates the performance of the method by applying the proposed controller to a battery emulator system and illustrates the results of the stability analysis. A discussion and a conclusion finalize this article. The optimality of the DMPC in unconstrained operation is shown analytically in the Appendix.

A. Overview
Several identical agents act in parallel on a central plant. The agents are not directly coupled; nevertheless, due to the feedback of the central plant, the agents affect each other indirectly. The main target is to control the central system output y c while compensating the effects of the measurable disturbance w, which has a direct influence on the central system. The central system has no direct control inputs and can only be manipulated by means of the agents (see Fig. 1).

B. Dynamic Model
The overall system consists of a set of M linear, dynamical actuators of order n a , and a linear central system of order n c . The number of inputs per actuator is denoted by m a . There is no direct state or input coupling between actuators, which is a necessary condition for the proposed DMPC method.
The dynamics of the i th actuator can be described by the linear, time-invariant state equatioṅ where A a ∈ R n a ×n a , B a ∈ R n a ×m a , E a ∈ R n a ×n c , and C a ∈ R 1×n a are the system, input, disturbance, and output matrix of the actuator, respectively. The disturbance matrix E a describes the feedback of the central system to the actuator. The actuator output z i is the i th input of the central system. The central system is described by the linear, time-invariant state equationẋ where A c ∈ R n c ×n c , B c ∈ R n c ×1 , E c ∈ R n c ×1 , and C c ∈ R 1×n c are the system, input, disturbance, and output matrix of the central system, respectively. The input matrix B c describes the influence of the actuator outputs on the central system. The set of actuators is defined by M = {1, . . . , M}.
The dynamics of the overall system, comprised of all actuators and the central system, is described bẏ where the system matrix A ∈ R (Mn a +n c )×(Mn a +n c ) , the input matrix B ∈ R (Mn a +n c )×(Mm a ) , the disturbance vector E ∈ R (Mn a +n c )×1 , the output vector C ∈ R 1×(Mn a +n c ) , the state vector x ∈ R (Mn a +n c )×1 , and the control vector u ∈ R (Mm a )×1 are composed as follows: ⎡ The coupling matrix G c = B c C a describes the effect of actuator states on central system states. The discretization of the continuous-time model of the overall system in (3) assuming a zero-order hold for the inputs and the sampling time T s gives the discrete-time model of the overall system

C. Load Distribution
The primary objective is to control the central system such that the output y follows the corresponding reference y ref , which generally varies over time. In addition, the disturbance w, which can be seen as load on the central system, must be compensated. In comparison to applying a single actuator, having several parallel actuators acting on the central system implies additional degrees of freedom. Restrictions are necessary in order to distribute the central system demands over the actuators in a controlled way. For this reason, references for the actuator outputs are defined. The reference z ref i for the i th actuator is based on the measured load by where κ i ≥ 0 are load distribution factors with i∈M κ i = 1. This ensures that the sum of all actuator outputs compensates for the load in steady state.

III. DECENTRALIZED MPC
In a DMPC scheme, the overall optimization problem is divided into several subproblems, which are solved in a decentralized manner. The control problem is optimized within the prediction horizon N p resulting in an optimal trajectory of control inputs. According to the receding horizon principle, only the first step of the optimal control trajectory is applied to the plant. At the next sampling instant, state measurements are updated and the procedure repeats.
Each actuator is equipped with a DMPC based on an appropriate subsystem model, which in combination is defined as an agent (see Fig 2). Since there are no direct agentto-agent couplings in continuous time, subsystems can be derived by merging the corresponding actuators with the central system. Provided that each subsystem has access to state measurements of all agents, the entirety of subsystems is an exact representation of the overall system in (3) ẋ However, the approach does not include the communication between agents. Therefore, the state information of other agents, i.e., the last term in (6a), cannot be considered. Not considering the influence of other agents entails problems. Each agent assumes to be solely responsible for the central system reference tracking and the load compensation. In fact, (M − 1) further agents contribute to the central system. This affects the control performance and provokes stability issues. Therefore, virtual subsystems are created, consisting of the corresponding actuator and a scaled version of the central systemẋ T , y i denotes the virtual central system output, andG c andẼ c denote scaled matrices.
Scaling includes two aspects: 1) the effect of the agent and the load on the central system is virtually amplified, which is equivalent to a reduction of the central system size, and 2) only a part of the load, which is defined by the agent's load distribution factor κ i , is considered. Consequently, each agent only deals with a virtual portion of the central system. Scaling principally targets the matrices G c and E c , which are amplified by M for the virtual subsystem,G c = MG c and E c = ME c . If the central system is a higher order system, additional scaling measures might be necessary. This includes postmeasurement scaling of certain central system states by κ i M. An example is given in the simulation results.
To be used by the DMPC, the virtual subsystem model in (7) is discretized by assuming a zero-order hold for the inputs and the sampling time T s , which gives The local objective function J i considers the deviation of the virtual central system output from the corresponding reference, the load distribution, and a penalty on the control inputs where ||·|| 2 Q denotes the square of the Q-weighted Euclidean norm. The matrices Q Y,D ∈ R N p ×N p and Q Z,D ∈ R N p ×N p are positive definite weighting matrices, and R D ∈ R m a N p ×m a N p is a positive semidefinite weighting matrix used for tuning. All agents apply the same weightings. The vectors Y i ∈ R N p ×1 , Z i ∈ R N p ×1 , and U i ∈ R m a N p ×1 denote the predicted trajectories of the virtual central system output and the i th agent output, and the input sequence of the i th actuator, respectively Note that predictions of the central system output Y i are virtual and vary between agents during transients if the load is unevenly distributed. The vector Y ref ∈ R N p ×1 is the output reference trajectory associated with (10), which is the same for all agents. According to (5), the agent output The i th agent minimizes J i in (9), which yields the optimal sequences of control inputs U i * where U i is the set defining the control input constraints of the i th actuator. The first step of the optimal control sequence is applied to the associated actuator. State constraints are not included because agent predictions deviate from actual trajectories in general.

IV. CENTRALIZED MPC
The CMPC works with the discrete-time model of the overall system according to (4). Measurements of all system states are transmitted to a central unit, which optimizes the global control problem and transmits the optimal control variables back to the actuators (see Fig. 3).
The objective function J of the CMPC has a structure analog to J i in (9) and is defined by where Q Y ∈ R N p ×N p and Q Z ∈ R M N p ×M N p are positive definite weighting matrices and R ∈ R Mm a N p ×Mm a N p is a positive semidefinite weighting matrix used for tuning. The set of DMPCs includes the output reference tracking term M times, whereas the CMPC does only once. In order to ensure J = i∈M J i , which is necessary for a fair performance comparison, the following relation between the reference tracking weighting matrices of DMPC and CMPC must be fulfilled: The matrices Q Z and R are block diagonal matrices with the main-diagonal blocks being Q Z,D and R D , respectively. This means that Q Z and R contain the same weightings as their DMPC counterparts. The vector Y ∈ R N p ×1 is the predicted output trajectory of the overall system. The contain sequences of the actuator outputs, the corresponding references, and the control inputs of all actuators, respectively Minimizing J according to (16) gives the optimal sequence of control inputs U * for all actuators where the set U covers the control inputs constraints of all actuators and O ∞ is the terminal set. The first step of each control sequence is applied to the corresponding actuator.

A. Optimality in Unconstrained Operation
The virtual subsystem model defined in (7) and the overall system model according to (3) are discretized assuming a zero-order hold for the inputs. Due to the scaling of the virtual subsystem, relations between blocks of the system matrix of the discrete-time virtual subsystem in (8a) (22) whereÃ is the system matrix of the continuous-time virtual subsystem, and blocks of the system matrix of the discrete-time overall system in (4a) can be found, namely, Similar relations can be shown for the control and disturbance input vectors. Note that mentioned relations between the subsystem and overall system matrices given here are only valid for the stated, simple central system dynamics. General extensions for higher order central system dynamics are out of scope here. An illustrative example is given in the simulation results.
Regarding the DMPC, predicted trajectories of the central system output Y i and the actuator output Z i can be expressed as where the matricesF,,,F Z , Z , and Z are derived from the discrete-time model as described in [26]. In contrast to W f , the disturbance trajectory W is shifted by one sample The matricesF = [F aFc ] andF Z = [F Z,aFZ,c ] can be decomposed into blocks related to actuator states and central system states.
The CMPC predictions of the central system output Y and the actuator outputs Z are defined analogously but include all actuators. Decomposing the corresponding matrices F, , , F Z , Z , and Z into actuator-related and central system-related blocks, similarly as in (23), one can show that aforementioned relations with respect to the subsystem matrices are maintained.
Assuming the unconstrained operation, the local DMPC objective functions and the global CMPC objective function can be optimized analytically. The optimal sequence of control inputs U * i for the DMPC is obtained by inserting (24) and (25) into the objective function J i according to (9) and setting the first derivative with respect to U i to zero, ∂ J i/∂ U i = 0. With (13), this results in Adding up all U * i over the agents considering Expressing the optimal sequence of all control inputs U * of the CMPC accordingly by analytically optimizing the global objective function J in (16) and considering the relations between subsystem and overall system matrices, such as F a = MF a , it can be shown that the CMPC solution for the sum of optimal control inputs i∈M U * i is exactly the same as for the DMPC in (28). Since all actuators are identical and linear, it can be shown that, consequentially, also the DMPC solutions of the resulting sum of actuator outputs i∈M Z * i and, therefore, all central system states, including the central system output Y * i , are exactly the same as for the CMPC. The complete derivation can be found in the Appendix.
Note that in general, individual agent inputs, states, and outputs differ from the corresponding CMPC solutions during transients and are therefore globally suboptimal. The DMPC achieves global optimality regarding individual actuator variables only in steady state because agents operate solely with measurements of the own actuator and the central system. In contrast, the CMPC possesses measurements of the overall system and operates globally optimal.
The load distribution can be adapted in steady state. In this case, the agent outputs converge to the adapted references preserving the optimality of the central system states. If the load distribution is changed during transients, the optimality is lost in general.

B. Stability Analysis
Assuming the output reference and the disturbance to be constant within the prediction horizon, i.e., . . , w] T , a linear state vector feedback control law for the optimal control input at the current step can be derived from (27) where the gainsk x ,k y , andk w are identical for all agents.
Combining the local control laws yields the linear state vector feedback control law of the DMPC in the centralized formulation Analogously, a linear state vector feedback control law can be derived for the CMPC Because of the central system scaling in the virtual subsystems of the DMPC, the gain vectorsK y and K y are identical. However, the feedback gain matricesK x and K x are different because the DMPC considers only the measurements of the own actuator and the central system for the computation of the control input of the i th actuator, whereas the CMPC considers the measurements of all actuators. Also, the gain vectorsK w and K w are different. Since the DMPC gains differ from the CMPC gains, individual control inputs resulting from (30) and (31) are different. Nevertheless, both feedback control laws yield the same sum of control inputs, as shown in Section V-A. Inserting the control law in (30) into the overall system model in (4) defines the closed-loop system dynamics of the unconstrained DMPC which includes an affine term as a function of y ref and w. The polytopic control input constraints U and physical state constraints X are defined by With the knowledge of the control law in (30), the control input constraints can be expressed in the form of state constraints, which allows to write a polytopic setP covering physical state constraints and control input constraints Note thatP depends on y ref and w.
HavingP, positive-invariant sets for the closed-loop system under the action of the DMPC can be derived. A setÕ ⊆P of initial states is called positive invariant for the closed-loop system (32) subject to the constraints in (35) if resulting trajectories will never violate those constraints at any time The maximal positive-invariant set or terminal setÕ ∞ ⊆P is defined to be the invariant set containing all positive-invariant sets inP [27]. This means that withinÕ ∞ , it is guaranteed that the DMPC operates without running into the control input constraints and without violating any state constraints. Consequently, the closed-loop stability of the DMPC is characterized by the eigenvalues of ( For asymptotic stability, all eigenvalues must have a magnitude strictly less than unity [28]. Regarding algorithms for computing the maximal positive-invariant set, the reader is referred to the literature [27], [29], [30].
Analogously, the terminal set O ∞ of the CMPC and the closed-loop stability of the CMPC within O ∞ are determined by the corresponding closed-loop system dynamics subject to constraints (33) and (34). Since the CMPC considers the measurements of all actuators and the central system in the optimization of the global objective function, the volume of O ∞ is expected to be larger than the volume ofÕ ∞ , which holds true in general. However, O ∞ is not a subset of O ∞ . In certain areas ofÕ ∞ , the CMPC runs into the control input constraints, whereas the DMPC operates unconstrained. The sum of control inputs of the unconstrained DMPC and CMPC would still be the same, but individual control inputs of the CMPC would exceed the boundaries.
Within the set R, which is the intersection of the terminal sets of the DMPC and the CMPC it is ensured that both the DMPC and the CMPC operate without violating any constraints. The stability of both control concepts is determined by the eigenvalues of the closed-loop systems. The sum of control inputs, the sum of agent outputs, and all central system states under the control of the DMPC are identical to the CMPC, as shown in Section V-A. Outside R, input constraints of the DMPC, the CMPC, or both potentially become active. Then, the sum of control inputs, and therefore the sum of agent outputs, and the central system states differ between the DMPC and CMPC in general.
To analyze the stability outside the terminal set, an N p -step admissible set can be computed for the CMPC, which has a global view of the system and includes terminal set constraints. The admissible set is the set containing all initial states that can be driven into the terminal set without violating constraints It can be developed in a backward step-by-step procedure beginning at the terminal set For more details, the reader is referred to the literature [27].
In contrast to the CMPC, the DMPC optimizes the control inputs locally and does not include terminal set constraints. Consequently, an admissible set cannot be derived for the DMPC. To analyze the stability outside the terminal set, the N p -step region of attractionÑ N p can be computed for the DMPC. In analogy to the admissible set of the CMPC, it is defined by the set of initial states that are driven into the terminal set in N p steps but under the optimal control inputs of the DMPC according to (15), which are optimized in each step The receding horizon principle must be considered, i.e., the whole sequence of optimal control inputs resulting from (15) cannot be applied directly, because of the deviation of the agents' virtual predictions from the actual trajectories. Based on (41),Ñ N p can be determined numerically by exploring the state space. A simple algorithm is proposed in Section VI-D.
Both the admissible set of the CMPC and the region of attraction of the DMPC assume the same constant y ref and w as in the computation of the corresponding terminal sets.

VI. SIMULATION STUDY
In this section, the application use case and controller setup are defined (Sections VI-A and VI-B), followed by a stability analysis of both unconstrained and constrained operation (Sections VI-C and VI-D). Then, a challenging test case is simulated and control performance is assessed both in unconstrained and constrained conditions (Sections VI-E-VI-G). Finally, test case extensions are discussed (Sections VI-H and VI-I).

A. Application
Consider the output stage of a battery emulator system for automotive testbeds. With the objective of increasing the emulator performance and handling higher load currents, two identical actuators are installed in parallel (M = 2) in order to control the central system (see Fig. 4). The central system consists of the capacitor C 2 , where the voltage v 2 is controlled by the actuators and the measured load current i p acts as a disturbance. The i th actuator consists of the filter capacitor C 1,i and a connection cable modeled in the form of a series connection of the inductance L 1,i and the resistance R 1,i . Each actuator equipped with a DMPC represents an agent. The demanded input currents i 1,i of the agents are provided by dc-dc step-down converter stages, one for each actuator. For the control of the converter stage, it is beneficial to use the time derivative of the actuator input current di 1,i/dt as control input instead of using i 1,i directly, and more details can be found in [31].
By defining the actuator state vector x i = [i 2,i , v c,i , i 1,i ] T , the actuator input u i = di 1,i/dt , and the actuator output z i = i 2,i , the linear state-space model of the i th actuator can be written aṡ where i ∈ {1, 2}. With the central system state x c = v 2 , which is also the central system output y c , and the disturbance w = i p , the first-order state-space model of the central system is given bẏ By defining the state vector x = [x 1 T , x 2 T , x c ] T and the control input vector u = [u 1 , u 2 ] T , the state-space model of the overall system can be written aṡ The system matrix in (44a) shows that the actuator states are not directly coupled. Nevertheless, the agents influence each other indirectly through x c .

B. Controller Setup
The subsystem model of the DMPC is derived according to (7). ConsideringG c = 2G c andẼ c = 2E c , the virtual subsystem model of the i th agent can be written as The weighting matrices are Q Y,D = 5 V −1 I, Q Z,D = 0.1 A −1 I, and R D = 0.1 ns A −1 I, where I denotes the identity matrix of appropriate order. The prediction horizon is chosen to be N p = 6. The CMPC, which is used as a benchmark to evaluate the performance of the DMPC, is based on the model defined in (44). The weighting matrix regarding the central system output tracking is scaled according to (17) giving Q Y = 10 V −1 I, whereas Q Z = 0.1 A −1 I and R = 0.1 ns A −1 I are only adapted with respect to the dimension.
The load distribution is chosen unevenly with κ 1 = 0.7 and κ 2 = 0.3 in order to demonstrate the capability of the DMPC.

C. Stability Analysis of Unconstrained Operation
By analytically optimizing the objective functions, linear state vector feedback control laws can be derived for the DMPC and the CMPC, see (30) and (31) The off-diagonal blocks of zeros inK x indicate that the DMPC does not consider states of other actuators for computing individual control inputs. One can verify that, apart from rounding errors in this rough presentation,K x x and K x x yield the same result if the states of both actuators in x are identical. This is a result of the central system scaling and shows that each DMPC implicitly assumes the states of other actuators to be the same as the own states for the consideration of coupling effects. With an uneven load distribution, the actuator states, however, differ from each other if the load is not zero. The arising deviation in steady state is compensated by the difference betweenK w and K w . Both control laws stabilize the system. As shown in Fig. 5, all eigenvalues of the closed-loop systems defined in (32) and (37) lie inside the unit circle of the complex plane.
Considering the control input constraints and physical state constraints, terminal sets can be computed for the closed-loop systems of the DMPC and CMPC for each pair of y ref and w. Figs. 8 and 9 show projections of the 7-D terminal sets onto the reduced 3-D state space of the two actuator outputs and the central system output for y ref = 100 V and w = 600 A. As the projections suggest, the volume of O ∞ is with 5.6 × 10 14 clearly larger than the volume ofÕ ∞ , which is 3.  For computing and plotting the polytopic terminal sets, the MATLAB-based Multi-Parametric Toolbox 3.0 was used [32].

D. Stability Analysis of Constrained Operation
Because of the high dimensionality of the system's state space, the step-wise procedure for determining the N p -step admissible set of the CMPC, see (40), is computationally extremely demanding and could not be executed. Therefore, N 6 was approximated numerically with a simple algorithm.
Step 1: Define an empty set V, where vertices of the admissible set will be stored. Define search directions and a starting point that satisfies the definition of the admissible set of the CMPC (39).
Step 2: From the starting point, go in a direction with a predefined step size until the first infeasible point is found, i.e., where (39) is not satisfied.
Step 3: Apply the bisection method over the range of the last feasible and the infeasible point to approximate the border of the admissible set. Add the last feasible result of the bisection method to V.
Step 4: Go to step 2 and repeat until all predefined search directions were explored.
Step 5: Compute the convex hull of the vertices in V to estimate the admissible set.
The algorithm was initially executed for the steady state defined by y ref and w as starting point and then repeated for all found vertices in a second iteration. The same algorithm was used to approximate the six-step region of attraction of the DMPC but using the corresponding definition according to (41).
Projections of the DMPC's region of attraction and the CMPC's admissible set are shown in Fig. 6. Both are supersets of the corresponding terminal sets. The CMPC's admissible set is clearly larger than the DMPC' region of attraction as expected. However, considering that the agents of the DMPC operate only with local measurements and do not include terminal set constraints, the DMPC's region of attraction covers a remarkable part of the CMPC's admissible set.

E. Simulation
The control performance of the DMPC is compared to the performance of the CMPC in a simulation covering four events.
1) The simulation is initialized with an output reference voltage of 20 V and a load current of 600 A but with the system not being in a steady state. While actuator 1 provides the target output current according to the specified load distribution with x init 1 = [420 A, 20 V, 420 A] T , the currents of actuator 2 are arbitrarily chosen too high with 2) The output reference voltage steps up from 20 to 100 V at t = 1.5 ms. 3) Between t = 3 ms and t = 4.2 ms, the load current changes dynamically with a gradient of ±500 A ms −1 . 4) At t = 5 ms, the load current increases abruptly from 400 to 600 A. Future trajectories of the output reference voltage and the load current are assumed to be unknown to analyze the interaction between agents and the closed-loop stability for unpredicted events. Current values are kept constant within the prediction horizon.
The results of the simulation are shown in Fig. 7. One can distinguish between unconstrained and constrained operations, the analysis of which is covered in the following. Note that the specified physical state constraints are not violated at any time of the simulation.

F. Performance Analysis of Unconstrained Operation
During events 1) and 3), both DMPC and CMPC operate unconstrained. Individual agent output currents of the DMPC differ slightly from the corresponding CMPC solutions during transients, as shown in the second plot of Fig. 7. Despite that, the sum of the agent output currents of the DMPC is the same as of the CMPC even though the local DMPCs do not include measurements of the other actuator's states. As shown in the third plot, the DMPC consequently achieves exactly the same control performance as the CMPC regarding the central system output (see Section V-A).
The fifth plot compares the objective function values of the DMPC and the CMPC throughout the simulation, where for the DMPC, the sum of the two local objective function values J 1 and J 2 is shown. At the beginning of transients, the objective function of the DMPC increases only to slightly higher values than its CMPC counterpart, indicating that the DMPC remains close to the global CMPC optimality. Afterward, the objective function values of both DMPC and CMPC decrease de facto exponentially over time, whereby the value of the DMPC converges toward the CMPC solution. In a steady state, where each agent compensates for the specified part of the load current, the DMPC operates with the global optimality of the CMPC.

G. Performance Analysis of Constrained Operation
Event 2), which is a large output voltage step, and event 4), which is a large load current step, represent two fundamental events in which the control input constraints get active.  The interaction between agents and the convergence during constrained operation are analyzed based on those two events. Fig. 8 shows the trajectory of event 2) in the reduced state space of the two actuator output currents and the central system voltage. The transition begins and ends on a plane determined by the specified steady-state load distribution The initial state is outside the terminal sets of DMPC and CMPC. With both controllers, the two actuators contribute equally to follow the output voltage reference, i.e., the deviations of the actuator output currents from the steady-state currents specified by the load distribution are the same. Consequently, the control inputs are set synchronously and therefore run into the constraints at the same time. The state-space trajectories of DMPC and CMPC coincide and lie within a plane perpendicular to the i 2,1 -i 2,2 plane. The plane is characterized by one of the complex eigenmodes of the closed-loop systems. With monotonically decreasing objective function values, see fifth plot of Fig. 7, both DMPC and CMPC drive the system into the corresponding terminal sets, where closed-loop stability is proven (see Section VI-C). The state-space trajectory of event 4) is shown in Fig. 9. Unlike the projected sets in the reduced state space indicate, the initial state lies outside the terminal sets of DMPC and CMPC. Initially, both controllers aim for compensating the higher load according to the specified load distribution, which results in different control inputs for the actuators. Since the load cannot be compensated instantaneously, the central system voltage v 2 drops. Consequently, the control task includes returning v 2 to the reference value besides compensating for the higher load current. Therefore, the resulting state-space trajectories of both DMPC and CMPC do not lie in K as they are affected by the aforementioned complex eigenmode.
The initial control input of actuator 1, which is higher than the control input of actuator 2 because κ 1 > κ 2 , runs into the constraint. The CMPC amplifies the contribution of actuator 2, see the first plot of Fig. 7, whereas the DMPC does not because agent 2 is not aware of the active control input constraint of agent 1. The sum of actuator output currents differs between DMPC and CMPC, and therefore, the DMPC does not preserve the CMPC optimality regarding the central system state, as shown in the third plot of Fig. 7. Nevertheless, both DMPC and CMPC decrease the objective function values monotonically and drive the system into the terminal sets, where the corresponding closed-loop systems are asymptotically stable (see Section VI-C).

H. Simulation With Predicted Events
The simulation results shown so far showed the closed-loop behavior for unpredicted events. If predictions of the output reference voltage and the load current are, however, available and utilized in the MPCs, the performances of both DMPC and CMPC improve significantly, in particular with respect to the output voltage. Fig. 10 shows the simulation results, including the same events, as in Section VI-E.
Because both controllers are able to prepare for the predicted events, control inputs stay within the specified constraints. The DMPC achieves the CMPC performance regarding the output voltage throughout the entire simulation, including the abrupt load step at t = 5 ms. Since the difference between actuator states of DMPC and CMPC is remarkably small, the overall DMPC performance is close to the global optimum, see the last plot of Fig. 10.

I. Higher Order Central System
In order to demonstrate the capability of the DMPC for higher order central systems, the central system of the application described in Section VI-A is extended by an LC filter with the parameters L 2 = 250 μH, R 2 = 0 , and C 3 = 4.6 mF (see Fig. 11). The two actuators used to control the central system are the same as in Section VI-A. The central system output is the voltage v 3 . Defining the central system state In addition, the intermediate current i 3 of the central system must be scaled after measurement,ĩ 3 = κ i Mi 3 . This is necessary to hold the balance between currents since the actuator output and the load current are virtually amplified through the central system scaling (see Section III). For closed-loop stability, the prediction horizon must be longer, N p = 15. Apart from that, the same controller settings and load distribution as in Section VI-B are chosen. The output reference and load trajectories are assumed to be predicted.
The simulation is initialized in steady state. At t = 3 ms, the output reference voltage steps up from 20 to 100 V, and starting from t = 10 ms, the load current increases from 200 to 600 A with a gradient of 500 A ms −1 . Despite the considerable deviation of the individual DMPC actuator output currents from the corresponding CMPC solutions, see the second plot of Fig. 12, the sum of the actuator output currents is the same for both controllers throughout the entire simulation. Consequently, the trajectories of the DMPC and the CMPC coincide regarding all three central system states, namely, v 2 , i 3 , and v 3 , as shown in the third and fourth plots. This demonstrates that the DMPC method achieves the global optimality of the CMPC regarding all central system states. In a steady state, the agents share the load current according to the specified load distribution.
Note that if the third-order central system is damped with R 2 > 0 , the entry on the main diagonal of A c with respect to the central system state i 3 must be scaled by (κ i M) −1 in order to derive the virtual subsystem model. The system matrix of Relations between subsystem and overall system, as discussed in Section V-A, are not applicable in this case. The optimality of the DMPC regarding central system states is lost. Nevertheless, simulations show that the central system performance of the DMPC still remains close to the optimum of the CMPC.
VII. DISCUSSION It has been shown that the set of DMPCs exactly achieves the optimal control performance of the CMPC with respect to all central system states. In the case of even load distribution, κ i = M −1 , all agents apply the same actuator control inputs. Since the sum of actuator control inputs is unchanged in comparison to the CMPC, it follows that in this case, the DMPC achieves the CMPC optimality also for individual actuator states and, thus, the overall system.
The incorporation of state constraints requires additional measures. As mentioned in Section III, virtual predictions vary between agents depending on κ i and deviate from actual trajectories resulting from the superposition of all agents' actions. The investigation of solutions is a potential topic of future research. If the load is evenly distributed and the agents operate in a synchronized manner, the consideration of state constraints is possible.

VIII. CONCLUSION
A DMPC method without communication for plants controlled via multiple, identical actuators is presented. Each DMPC is based on a virtual subsystem model consisting of the corresponding actuator and a scaled version of the central system. In this way, each agent, an actuator equipped with a DMPC, virtually deals only with a portion of the central system. An analytical derivation shows that in the unconstrained operation, the set of DMPCs generates exactly the same sum of control inputs as an equivalently tuned CMPC and that the DMPC consequently achieves the optimal performance of the CMPC with respect to all central system states. Agent outputs generally differ from corresponding CMPC solutions during transients but converge afterward. If the load of the plant is evenly distributed over the agents, the DMPC operates with the overall optimality of the CMPC. Stability in unconstrained operation is determined analytically.
In constrained operation regarding the control inputs, the stability is analyzed by computing the region of attraction, which is compared to the CMPC's admissible set. The interaction between agents and the convergence is investigated based on numerical studies of a battery emulator system. The simulation results of this application show that the DMPC achieves the central system optimality for a range of events in ideal operation, where predictions are accurate. If abrupt transients are, however, not predicted, the DMPC preserves stability and operates closely to the optimum considering the control input constraints.

APPENDIX A OPTIMALITY OF THE SUM OF CONTROL INPUTS IN
UNCONSTRAINED OPERATION The local objective function of the i th DMPC [see (9)] is (13)]. The predicted trajectories of the central system output Y i and the agent output Z i can be written as where the matricesF,,,F Z , Z , and Z are derived from the discrete-time model as described in [26]. The matrices F andF Z are shown decomposed into blocks related to the actuator and the central system states. Inserting the expressions for Y i and Z i and analytically optimizing the local objective and finally results in the optimal sequence of control inputs U * i for the i th agent Adding up over all M agents considering i∈M κ i = 1 yields the sum of optimal control inputs of the DMPC The objective function of the CMPC [see (16)] is where the weighting matrices Q Z and R can be expressed by means of their DMPC counterparts and Z ref can be written as follows with K i = κ i I: The predicted trajectories of the central system output Y and the actuator outputs Z can be expressed analogously as for the DMPC but including all actuators . .
. . Adding up all rows on both sides by left multiplication of I I · · · I and considering the relation between blocks of the CMPC matrix Z and the DMPC matrix Z , namely, Considering further relations between CMPC and DMPC matrices, namely, F Z,a + (M − 1)F Z,o =F Z,a , a = 1 /M, F a = 1 /MF a , F c =F c , F Z,c =F Z,c , = 1 /M, and Z,a = 1 /M Z ; moreover, i∈M K i = I and (17), The sum of optimal control inputs of the CMPC is finally which is exactly the same solution as for the DMPC.

APPENDIX B OPTIMALITY OF THE SUM OF AGENT OUTPUTS IN UNCONSTRAINED OPERATION
The trajectory of agent outputs resulting from optimal control inputs is The trajectory of the sum of optimal agent outputs only depends on the current states, the disturbance trajectory, and the sequence of the sum of optimal control inputs, but not on the sequence of individual control inputs. As i∈M U * i is the same for DMPC and CMPC, it follows that also, i∈M Z * i is the same for DMPC and CMPC.

APPENDIX C OPTIMALITY OF THE CENTRAL SYSTEM STATES IN UNCONSTRAINED OPERATION
The trajectory of the central system output resulting from optimal control inputs can be written as